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ABSTRACT 


A numerical code, USPOTF2, has been formulated to solve for the potential flow 
for two airfoils executing unsteady motions in an inviscid incompressible flow mediuin. 
This code is an extension of an existing code U2DIIF, which does the same calculations 
for the single airfoil case. The technique uses the well known Panel Methods for steady 
flow and extends it to unsteady flow by introducing a wake model which creates a non- 
linear problem due to the continuous shedding of vortices into the trailing wake. The 
presence of the second airfoil introduces a set of non-linear coupled equations for the 
Kutta condition. Numerous case-runs are presented to illustrate the capability of the 
code. The case of the step change in angle of attack is compared with Giesing’s work. 


All other case-runs are illustrated together with the results for the single airfoil case. 


THESIS DISCLAIMER 


The reader is cautioned that computer programs developed in this research mav not 
have been exercised for all cases of interest. While every effort has been made, within 
the time available. to ensure that the programs are free of computational and logic er- 
rors, they cannot be considered validated. Any application of these programs without 


additional verification is at the risk of the user. 
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I. INTRODUCTION 
A. GENERAL 


In this paper , a numerical method is formulated to solve for the flow about two 
two-dimensional airfoils which are arbitrary located and are performing an arbitrary time 
dependent motion in an inviscid incompressible fluid. The original work by Teng [Ref. 
1] is for a single arbitrarily defined airfoil in the same potential flow condition. The ex- 
tension of Teng’s code to two bodies is considered here. Where possible, the same con- 
vention and notation are adopted. As for the single airfoil case, all velocities are 
non-dimensionalised with respect to the free stream and all lengths with respect to the 
chord length. 

A new subroutine (SUBROUTINE NEWPOS) is added to transform either of the 
two local coordinate svstems to the global coordinate svstem. A modification 1s entered 
into the treatment of the Kutta condition to make it consistent with the unsteady flow 
Kutta condition treatment. Subroutine COFISH 1s deleted and its role is taken over bv 
Subroutune COEF which normally performs the formulation for the flow tangency con- 
dition for the unsteady case. A more accurate method 1s also introduced to obtain the 
Select wmaotential. lhe reader is referred tothe work of ‘Krainer (Ref. 2] for further 
improvement of the original code. 

This documentation is set up as for the original documentation 1n order that it will 
be easier to follow and cross-referenced. While it is the intent of the author to keep this 
thesis as complete as possible. the reader is well advised to review Reference | for the 
work involved in the single airfoil case: no special effort will be made to reproduce it 


here. 


B. BASIC THEORY AND APPROACH 
1. Steady Flow Problem 
The treatment for the two airfoils case in steady incompressible flow follows 
closely to the single airfoil case. The governing equation, as for the single airfoil case , 
follows from the Conservation of Mass and the Condition of Irrotational Flow. 
The continuity equation of an incompressible fluid (div V= 0) and the condi- 


tion of irrotational flow (curl V= 0) leads to the well known Laplace Equation. 


div( gradd ) = 0 (1.1) 


with @ denoting the disturbance potential for the velocity. This is seen to be the classic 
Neumann problem of potential theory with the usual problem of defining the boundary 
conditions. The boundary conditions for the disturbance potential @ are that its gradi- 
ent normal to the surface be equal to the normal velocity of the surface and that its 
gradient vanish at infinity; that is 


Vien =V.enonS (1.2) 
lim V4( P) + 0 (1.3) 


where ie is the resultant velocity of a point on the body as seen from the inertia frame 
of reference, 7 is the outward unit vector normal to the body, S is the bodv surface and 
P represents a general point. Equation 1.2 holds for both airfoils i.e. on both surfaces. 

The pressure coefficient is obtained through the Bernoulli’s equation which de- 
rives from the Momentum equation. 

The approach adopted 1s associated with Hess and Smith [Ref. 3] who devised 
the popularly known PANEL method in the early sixties. In words, the boundary or 
airfoil surfaces S, and S, about which the flow is to be computed is approximated by a 
large number of surface elements whose characteristic dimensions are small compared 
to those of the body. Over each surface element, a uniform source distribution and a 
uniform vorticity distribution is placed. The source strength (q,) varies from element to 
element, while the vortex strength (7,) is the same for all elements in the same airfoil but 
is different across the airfoil. The singularity strengths are determined from the flow 
tangency condition on both body surfaces and the two Kutta conditions at both trailing 
edges. With the determination of the singularity strengths, the relevant aerodynamic 
data can then be subsequently computed. 

2. Unsteady Flow Problem 

The unsteady problem is similar to the steady flow problem in that thev both 
have the same governing equation viz. the Laplace equation, and that for both problems, 
the pressure and velocity are decoupled so that the velocity and pressure calculation can 
be computed separately and consecutively. 

This problem differs from the steady flow in that another model is required to 
simulate the continuous shedding of vorticity into the trailing wake. The existence of a 
vortex sheet behind the airfoil can be explained by the Helmholtz theorem which 1s 


basically a statement of the Conservation of Vorticity. This requires that anv change in 


the circulation around the airfoil must be matched by an equal and opposite vortex 
somewhere in the flowfield. The presence of the countervortices provides the flow with 
a kind of a memory in that the flow at a particular time is affected by the bound circu- 
lation of the past. It is this non-linearity that distinguishes the numerical technique from 
the simple steady flow problem of solving N linear equations in N unknowns. 

The solution technique requires an iterative type solution. The present ap- 
proach follows closely the original panel method of Hess and Smith as described in the 
steady flow development, while with regard to the modelling of the wake, it adopts the 
procedure advocated by Basu and Hancock (Ref. 4]. A uniform source distribution (q,), 
and a uniform vorticity distribution [7(4], as for steady flow is placed on each panel at 
time s, where } denotes the panel number and / the airfoil number. The wake consists 
of a single vorticity panel attached as an additional element on each airfoil through 
Which the vorticities are shed into the respective wake and a series of point vortices 
which are being convected downstream with the fluid. A uniform vorticity distribution 
of strength (:,(/)), 1s placed on the wake panel of each airfoil. This panel is further 
characterised bv its length A(/, and its inclination ( ©,_),), With respect to the respective 
local frame of reference. After each time step. the vorticity of the wake panel 1s con- 
centrated into a single point vortex and convected downstream. Simultaneously a new 
wake panel is formed. The downstream wake of point vortices is thus formed bv the 


shed vorticity of previous time steps. 


C. SCOPE 

Chapter I] extends the original code to handle the steady flow problem for two 
airfoils set at different relative distances and angles of attack. 

Chapter III deals with the unsteady problem for the two airfoils svstem. It intro- 
duces a new subroutine and a more accurate method of calculating the velocity poten- 
tial. The Kutta condition for the two airfoils system 1s specially treated as its unique 
problem of a non-linear coupled system is not seen in the single airfoil case. 

Chapter IV describes the computer program, its essential capabiliues and limita- 
tions, its associated subroutines, its input requirements and its associated output print- 
Out. 

Chapter V presents the results of some case-runs. Of interest is a Comparison case 
of a step change in angle of attack (AOA) with Giesing [Ref. 5] for the same airfoil 
undergoing an impulsive start at the same AOA. Case-runs will also be run with both 


airfoils at large distances apart to compare with the single airfoil case. In addition, ex- 


ample cases for which no comparisons exist are given, to exhibit the capabilitv of the 
method. 
Finally, Chapter VI concludes with future development efforts and the application 


potential of this numerical method. 


Il. STEADY FLOW PROBLEM FORMULATION FOR TWO AIRFOILS 
A. GENERAL 


The modification work for the steady flow is straightforward. The revised program 
allows for two arbitrarily defined airfoils placed at an arbitrary distance set at different 
angles of attack. For reason of simplicity, the number of panels and nodes and the pivot 


location are set to be the samie for both airfoils. 


B. FRAMES OF REFERENCE 

Three frames of reference are involved in the two two-dimensional airfoils’ case in 
steadv flow. These are three inertia frames of references as indicated in Figure 1. 

The first inertia frame of reference (also known as global frame of reference) is set 
at the pivot position of the first airfoil with the X-axis pointing in the direction of the 
free-stream velocity. The two other inertia frames of references. henceforth, will be 
known as two frozen local frames of reference! (x,1, and x.1;) are fixed respectively to 
each airfoil with the x-axis coinciding with the chord line originating from the respective 
leading edge. The two local frames of reference are set apart by AShift and YShift on 
the global frame of reference. In steady flow, the fluid velocity and pressure depend onlv 
on the spatial coordinates (X.Y) and not on time. 

The two airfoils are defined in the local coordinate svstem as input data for sim- 
plicitv and are then transformed to the global coordinate svstem through a knowledge 
of the relative positions of the 2 airfoils’ pivots positions and the respective local angles 


of attack. 


C. STEADY FLOW PANEL METHODS 
1. Definition of nodes and panels 
Each airfoil surface is divided into n straight line segments called panels bv 
(n+1) arbitrarily chosen points called nodes. The numbering sequence begins with 
panel 1 on the lower surface at the first airfoil trailing edge and proceeds clockwise 
around the airfoil contour so that the last panel on airfoil 1 ends on the upper surface 
at the trailing edge. This numbering sequence then proceeds in a similar fashion for the 


second airfoil ( See Figure 2). As with the single airfoil case, the numbering sequence 


1 For the steady case, the notation ‘frozen’ will be dropped 


a XSHIFT | 
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Figure 1. Frames of Reference for Steady Flow. 
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dictates that the airfoil body always lies on the right hand side of the i panel as one 
proceeds from the 7” node to the (¢+ 1)" Also the 1" and (m + 1)" nodes coincide at the 
first airfoil trailing edge and the (+ 2)" and (27 + 2)" coincide at the second airfoil 
trailing edge. This numbering system facilitates the definitions of the unit normal vector 
n, and the unit tangential vector /, for all panels with n, being directed outward from the 
body into the flow and z, directed from the /* node to the (i+ 1)" node. The numbering 
of the panel svstem is somewhat complicated by the fact that a continuous panel num- 
bering sequence across the airfoil is desired. This procedure leads to the peculiar panel 
numbering behaviour seen in Figure 2 where for the first airfoil, the i panel lies between 
ieoeeeti + 1)" modes while for the second airfoil, the i* panel lies between 
(7+ 1)" and (7 + 2)* nodes?. 
2. Distribution of Singularities 

Over each surface element of the two airfoils, a uniform source distribition and 
a uniform vorticity distribution is placed. The source strength gq, varies from element to 
element for each airfoil while the vorticity strength y, remains the same for all elements 
in the same airfoil but is different across the airfoil. This choice of singularities follows 
closely the original panel method of Hess and Smuth. [t automatically satisfies the 
Laplace Equation (which is the governing equation for the inviscid incompressible flow) 
and the boundarv condition at the far field (co). In addition, as the Laplace Equation 
is a linear homogeneous second order partial differential equation, an overall compli- 
cated flow field can be built up by the combination of simple flows with the condition 
that the appropriate boundary condituon on the airfoil be satisfied accurately. 

For our case, the overall flow field (represented by the velocity potential @) can 
be built up by three simple flows. Writing this in terms of the respective local frame of 


Revererice. 
D(x 3) = O(% 1H) + Os(x,¥) + OX ,¥) (2.1) 
where @ (x , 4) 1s the potential of the onset flow, 


}..(x .y¥) = VU [xCos a() + ySin o(/] (222) 


2 For the code 6, has been defined out of the above convention with 6,,i1= 1,2..n as per panel 
numbers for the first airfoil. but with @,_, reserved for the wake element of the first airfoil and @,, 
i=n+2n+3...2n+] for the second airfoil with @,,_, reserved for the wake element of the second 
airfou for unsteady flow. 
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Panel j : 
Source distribution = q; 
Vorticity distribution = 7; 
Notie: 1.Nodal points! defined by i = 1320). ntl Sfom first sage 
and i = n#+2 ... 2n+2 for second airfoil. 
2.Panel number defined by j = 1,2... n for first airfoil 
and j = n+l ... Zn forgsecongeair1ou, 


Figure 2. Panel Nethods Representation fur Steady Flow. 


@, is the velocity potential of a source of distribution q(s) per unit length. 


le — feuds (2.3) 





?, 1s the velocity potential of a vorticity distribution y (s) per unit length 


o,=- |= 6 ds (2.4) 


At this point, the disturbance potential, @. is introduced, which is defined to be the sum 


of the potential due to the source and vorticity distribution. 

9=9¢,+9, (2.5) 
Equation (2.1) can then be read as 

b=o +6 (2.6) 


The convenience of defining the above allows for the total velocity vector to be viewed 
as two components viz. the onset and the induced velocity due to the disturbance po- 


tential. The total velocity is thus: 


It is in the introduction of the disturbance potential that leads us to the concept of in- 
fluence coefficient which will be elaborated in a later section. 
The pressure coefficient can be obtained from Bernoullt’s Equation which 1s 


derived from the Conservation of Momentum. 


<2 = y — J eset (2.8) 
— pVe l -_ j 


3. Boundary Condition 
As in the single airfoil case, the boundary conditions to be satisfied include the 
flow tangency conditions and the Kutta condition. The flow tangency conditions are 
Satisfied at the exterior mid-points (control points). The normal velocity is taken with 
respect to the respective local frame of reference (for consistency with unsteady flow 


notation to be introduced later): 
();=0, i= 22ers eS) 


where (V*), is the normal component of the total velocity3. 

The Kutta condition postulates that the pressure on the upper and lower sur- 
face at the trailing edge of each airfoil be equal. For steady potential flow, equal pres- 
sure implies equal tangential velocity in the downstream direction at the first and last 
panel of each airfoil viz the Bernoulli equation. With our definition of the tangential 


vector we then have 
(i), .-(V), lst airfoil 
(HY), = — (Vy, Quid airfoil (2.10) 


As with the single airfoil, equations 2.9, and 2.10 lead to a linear svstem of 
(2n+ 2) simultaneous equations. With "and V" expressed explicitly in terms of ¢, (j = 
1,2,..n. n+l... 2n) andy,(l=1, 2) Wehave 2n + 2Zunknowns ind 2n + 2)isisiem Gilinean 


Simultaneous equations which can be easily solved. 


D. INFLUENCE COEFFICIENT 
1. Concept of Influence Coefficient 

As introduced earlier, this important concept of influence coefficient results 
from the presence of the disturbance potential which follows from the presence of the 
singularities. Formally, an influence coefficient 1s defined as the velocity induced at a 
field point by a unit strength singularitv (be it a point singularity or a distributed 
singularity) placed anywhere within the flow field. Recall that equations 2.9 and 2.10 
require the computation of the normal and tangential velocity components at all the el- 


ements control points. The normal components of velocities are essential in satisfying 


3 This follows Giesing’s notation where the velocity with respect to the airfoil frame is denoted 
as total velocity. While there is no misunderstanding in the steady flow, the notation actually refers 
to the relative velocity with respect to the airfoil moving frame of reference for the unsteady flow. 
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the flow tangency conditions, while the tangential components of velocities are necessarv 
for satisfving the Kutta condition as well as computing the pressure distribution. 
2. Notation for Influence Coefficient 

The notation used in this documentation will be as for the notation for the 
single airfoil case with a slight modification to take into account the effects of the the 
second airfoil. As the influence coefficients are related to the geometry of the airfoil and 
their relative positions, it must, of necessity be computed with respect to the global 
frame of reference for the two airfoils’ case. 


For steady flow, the following influence coefficients are defined. 


A? : normal velocity component induced at the i” control point by unit strength 
source distribution on the /™ panel. 


e A! : tangential velocitv component induced at the # control point by unit strength 
source distribution on the /” panel. 


e 8B: normal velocity component induced at the /“ control point by unit strength 
Vorticity distribution on the j” panel. 


e 8B : tangential velocity component induced at the 7” control point by unit strength 
vorucity distribution on the j” panel. 


where 1 and j denotes panel numbers and are defined as : 
mee ie. ..n, atl... 2n 
ae), ...n, n+! ... M 
3. Computation of Influence Coefficient 


A single source located on the j“ panel in the local frame of reference (see Figure 


3) induces a total velocity of 





aS ue 


where the component perpendicular to the inducing panel 1s 


h 








A an ae —s a y | 

D2 ie a Cale, 

and the component parallel to the inducing panel 1s 

yi. (-+)= eS (2.13) 
Zar r De Aen 


For a distributed source panel on the j” panel in the local frame of reference, we have 


It 


=a Sy 2 J 


(jar > Yar) 





Figure 3. Influence Coefficient for Source Panel 


53 

, ‘ik l 

V5 = |. V(s)ds = 5 FIAN (2.14) 
! 


bo 4 





where 
FTAN = tan7 (xem, — x;)Qvity = Yj4.4) — Gory — YI; — X41) 
(xin; — x;)(xM; — X44) + QI — YM: — Yar) 
xm, = 1/2 (x; + %,,) 
ym, = 1/20; +41) ele) 
5 = 
iss “|. V,(s)ds = —>— FLOG (2.16) 
where 
cit 2 i ea meas + OS yen) (2.17) 
Cag Ge Gury 4) 
Transforming the above to global frame of reference, we have 
Ai; = V77Cos(8; — @;) — ¥7,Sin(9; — 4) (2.18) 
Ai, = ViSin(@,— 0) + V,Cos(8; — 4) (2.19) 
B; = — (17;Cos(6; — Oi) V’,Sin(@, — @))) (2720) 
Bi; = V,Sin(6,— 6) — V;,Cos(@; — 8) 21) 
For the Steady Flow code, we define 
AAN(I,J) = Aj = By (2.22) 
BBN(I.J) = — Ay = B; (223 
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E. NUMERICAL SOLUTION SCHEME 
1. Rewriting the Boundary Condition 
Though it is not critical for the steady flow problem to adopt any particular 
frame of reference, we will adopt the local frame of reference in satisfving the boundarv 
condition for consistency with the treatment used in the unsteady case. 

Even though the influence coefficients are defined in terms of the global frame 
of reference, the computation of the normal and tangential velocities of each panel is 
independent of the coordinate system used. Thus the normal and tangential velocities 
obtained with the global coordinates system would be the same as those obtained for the 
local coordinate system. Using the local coordinate svstem, the flow tangency condition 


of equation 2.9 is defined as follows : 


eS) = y, Aiq; + y fine HOY + V_ Sinfo() -—@)] = 90 (2.24) 


j=n+! j=n+) 


where 


al te oe 
— il on 


The Kutta condition of eqn 2.10, in terms of influence coefficients, becomes for airfoil 
ie 


¥, Avg, t oy Ayiq; + oe ay a+r, Cos[a(1) — 6,] = 


j=nt+l j=ntl 
2n 
-{) 4) oF 4 Anya + v1 oy By a y(2) » Bry a5 Us Cos[a(1) = sah (2:25) 
j=n+! j=n+!) 


and for airfoil 2: 


¥, Anat jj + ay Anti jG + Dy, Bens y(2) Saat Ves Cos[o{2 Br) = 


j=ntl j=nt+) 


2n 


yer Sy Aon jj SB ‘Pius se) » Bony ll eosin, |. (2.26) 
i=] 


j=n+l jJ=n+!) 


2. Solving for the Strengths of Source and Vorticity Distribution 
Equations 2.24, 2.25 and 2.26 can be written as a set of 2n+2 linear simultane- 
ous equations with 2n+2 unknowns (q,,j= 1,2... 2” and y,,/= 1,2) and solved as for 
the single airfoil case. However, to make the routine consistent with the unsteady flow 
case, the following method 1s adopted. 
The flow tangency condition can be rewritten explicitly for the g, in terms of 


y(1). (2) and the free stream constant term; that is 


peep? Cpa oe * ¢ | Gyn q} Q 2n4) Q) 2n+2 by 
Cay Sie ee ° ° *). 22,2n q2 Q2 2n4+1 Q 2n42 by 
e322 |) 330% * * «32» q3 Q3 2n4+1 Q3 2n+2 b; 
3 P oe ee a — ; y(1) + ; “(2)+] « 
Steere fins ° * ° -82n2n |) Fn @2n,2n+! @2n,2n+2 ban 


Gauss Elimination is then used to solve for the qg, in terms of y(1), }(2) and the 


constant term. This gives 


q, = bi,y(1) + 62,y(2) + 63, f = 1,2... 20 (2.28) 


Equation 2.28 is then substituted into the Kutta condition at the two trailing edges to 
form two linear simultaneous equation with two unknowns y(1)andy(2) . Gauss 
Eimination is again used to solve for the Vorticity distribution y(1) and y(2) and these 
are then back substituted to solve for the q,. 
3. Computation of Velocity and Pressure Distribution 
Once the q, (j = 1.2... 2) and y, (/= 1,2) are solved, the velocities at all the panel 


control points can be easily obtained. The normal velocity is given by 


is 


n 2n 
(V"), = » Ajay yy Ai}4j +y DY: By + 12))- By + Fo. Sin{ a() — 9,] (2729) 
j=! 


j=n+l j=nt! 


while the tangential velocity is given by 


n 2n n 2n 
(V'), = » Aig a0 » Aig; + (1) By + y(2) » By + V,, Cos «(2 — 4] (2.30) 
J! j=ntl je j=nt] 


with 
f=l2 nt lle 


= ole oh ae 


The normal velocity will obviously satisfy equation 2.9 (used to check the code) showing 
that the flow tangencv condition 1s satisfied while the the total velocity will be given by 


the tangential velocity. 
, Bi : 
} fala (J Nips i= L.2@ 2n (2.31) 


where 


27 


--)y Big, Sy Brg, + HOLA (2) ) AP + V,,Cosl a()—6;) (2.32) 


j=ntl j=n+l 
Substituting equation 2.31 into 2.8 for C, with (V'), defined as in equation 2.32, the 
pressure coefficient at the /“ control point is 


(C,), = 1-(¥ fp, i= 1,2...20 (2.33) 


4. Computation of Forces and Moments 
The two-dimensional aerodynamic lift (C,), drag (C,) and pitching moment (C,) 
are calculated with respect to the global frame of reference. The moment coefficients 


are computed with respect to the respective leading edges. For the code, we have 


Ck) = SANE. — X}) (2.34) 


an = Sic ro. — 1) 


Cull) = > (Cpl (Sigs — XDA + (Yigg — YOY ny) 


where / = 1,2 and X, X_,, ¥,, ¥..,. Am, Ym, are defined as before. 


Wi. UNSTEADY FLOW PROBLEM FORMULATION 


A. GENERAL 
The modification work for the unsteady flow is much more involved. It requires the 
following: 

1. The establishment of five frames of reference viz one fixed inertia frame of reference 


(global), two moving local frames of reference4 and two frozen local frames of ref- 
erence.> 


Reformulation of the two Kutta conditions which are coupled non-linearly. The 
solution requires an iterative procedure to compute for the two y(/). There are two 
possible solutions to the Kutta condition due to the quadratic nature of the 
equations. The solution which ensures that the product of the tangential velocities 
of the first and last panels of each airfoil 1s negative 1s accepted as the solution. 


. The creation of a new subroutine (SUBROUTINE NEWPOS) which transforms 


all coordinates in either of the two respective local frames of reference to the global 
frame of reference. This simplifies the definition of the airfoil, wake element and 
core vortices relative global geometries. The code requires that the airfoil be de- 
fined with respect to the respective frozen local frames of reference once vomee 
Subsequent time dependent airfoil motion, wake panel behaviour and core vortices 
convection are computer generated. 


. The introduction of a more accurate method to obtain the velocity potential bv 


integrating the velocity over smaller panels on the airfoil without having to store 
large arravs of influence coefficients which are not needed for satisfying the flow 
tangencyv condition. 


Extension of the influence coefficient to include the effects of the second airfoil with 
its oWn peculiar weke. This also requires an introduction of an additional influence 
coefficient, that on the wake element due to the wake element from the other 
airfoil. 


B. FRAMES OF REFERENCE 


The inertia (global) and the two frozen local frames of reference at a specified time 


1, are defined as for the steady flow case. The two moving local frames of reference have 


their x-¥ axes as for the frozen local frame of reference but this frame 1s moving with the 


airfoil. 


4 Moving local frames of reference are used to satisfy the flow tangency equation which sim- 


plifies to equation 2.9. 


5 Frozen local frames of reference are inertia frames (used in the steady flow case) of reference 


used to convect the core vortices of the previous time steps with respect to that time step frame of 
reference and subsequently transformed to the current time step frozen frame of reference. 
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C. UNSTEADY FLOW MODEL 
1. Rigid Body Motion 

The rigid body motion for the two airfoil svstem is an extension to the single 
airfoil svstem. Both airfoils are considered to have a mean velocity of — V’_, with a time 
dependent translational velocity of — [UD,i + Vj] and a rotational velocity — Q(), 
which can be in phase or out of phase; i and j are unit vectors in the respective local 
frames of reference and 22 is positive in the clockwise direction. The flow will be deter- 
mined with respect to the moving local frame of reference. The flow tangency condition 
takes its simplest form in the moving frame of reference. The flow tangency condition 
seen from this frame of reference will satisfy equation 2.9. The unsteady stream velocity, 
Vimeem 1S Made up by the vector sum of a mean velocity V,, a time dependent 


translational velocity (U(),i - (0 ] and a rotational velocity Q(/, (Figure 4). 


Veoream = Vz,{ Cos| o(D, Ji + Sinl (Dy i} + [UOg + MOG 1+ 20,67 -) BN) 


The disturbance potential, of necessity, is defined with respect to the inertia frame of 
reference, for it is onlv in this frame that the flow is irrotational. It is then transformed 
to the moving frame of reference for ease in treating the flow tangencv condition. The 
disturbance potential is redefined to include the contributions from the wake panel and 


the core vortices from both airfouls. 
6= 6,4 6,4 byt bo Ce 


Wie vera) velocity 1s then written : 


— a 


| etat = OF eter Bs Vo (52) 


The unsteady Bernoulli's equation for the pressure coefficients on the airfoil surface is 
written with respect to the moving frame of reference. Giesing showed this to be written, 


in our notation as: 


aa IY cece 2 aa 2 2 CH 
ee eas gata ye 3.4 
C, loge ( V a V ) ye Cl (3.4) 
Ds a 
where V,,,,, and V,,,,, are defined according to equations 3.] and 3.3 respectively. 


iy 





xe 
e) 
| 


Vitdy 


V stream = Vyo{ Cos a(l)yi + Sin o(1)yi} + (UD gi + Vi } + QW), 07 — xf) 


Figure 4. Frames of Reference for Airfoil 1 
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2. Wake Model 
Recall that the unsteady flow model requires an additional model to simulate 
the continuous shedding of vorticity into the trailing wake. The treatment of this vortex 
shedding process follows the approach of Basu and Hancock. The shed vorticity takes 
place through a small straight line wake element attached as an additional panel to the 
trailing edge of each airfoil with a uniform vorticity distribution [y,()],. This shed 
Vorticity panel will be established if its length A(/, and inclination O(/), to the respective 


local frames of reference satisfv the Helmholtz theorem : 
Are + 0 = yD (3.5) 
or AD clr le =F ei — PO = SSL Opa — 204] (3.6) 


Where SS 1s the perimeter of the airfoil and T,_, and y,_, are respectively the total circu- 
lation and vorticity strength which are determined at the previous time step J,_,. 

At the next time step 4,_,, the shed vorticity panel will be detached from the 
trailing edge and will convect downstream as a concentrated free vortex with circulation 
A(!), [7.(9). or F,.(9 — (0 at the resultant local velocity of the fluid particle. This wake 
convection process 1s illustrated in Figure 5 where the airfoils’ subscripts are dropped 
without loss of generality. 


In the code, the convection of the core vortices is broken into three steps: 


1. The core vortices are first convected using the resultant absolute velocity with re- 
spect to the frozen local frame of reference of the previous time step. 


ty 


This is followed by a transformation to the current frozen local frame of reference. 


3. Finally this is transformed to the global frame of reference by using the new sub- 

routine (SUBROUTINE NEWPOS). 

3. Additional Boundary Conditions 

The unsteady flow model has now introduced an additional boundary condition 
viz. the conservation of vorticity (equation 3.5 or 3.6) through the modeling of the wake. 
However, the introduction of the wake creates three additional unknowns for each 
airfoil, that is, y,(9,, A(D, and O(),. As such two additional conditions are required for 
each airfoil in order to solve the svstem. The approach suggested by Basu and Hancock 


is extended to the two airfoil case. 


1. The wake panel is oriented in the direction of the local resultant velocity at the 
panel midpoint. 
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Figure 5. 


Vortex Shedding at Time Step t, 


Helmholtz’s theorem 


Oy, (YW + T= Ty. 
Panel j 
Source Distribution (93), 


Vorucity Distribution Vx 


k 
(V..) 
Ven tanO, = —+& 
( Je k oh 
0, DF a a) aie 


(ep 


Panel Methods Representation for Unsteady Flow 


| 


Ou 


we (1), 
peo n= on 3.7) 


Where [}2()], and [F3()], are the x and y velocity components at the midpoint of 
the wake panel with respect to the frozen local frame of reference®. 


tY 


The length of the wake panel is propotional to the magnitude of the local resultant 
velocity at the panel midpoint and the step size of the time step. 


Ay = (te — tea IRON) + (VO? (3.8) 
We now have all the necessary equations to solve for the system. 


D. INFLUENCE COEFFICIENTS 
As with the boundary condition, additional influence coefficients need to be defined 
as a result of the wake model. The definitions in the single airfoil case are thus extended 
as follows: 
1. Mfore A’s and B’s Influence Coefficients 
Defining NP3 and NP-4 as the panel number for the wake for the first and sec- 
ond airfoil respectively and h as a arbitrary core vortex, the following additional influ- 


ence coefficient are required. 


e (B.,;), : normal velocitv component induced at the / panel control point by unit 
strength vorticity distribution on the wake panel of the first airfoil at time 4,. 


e (B:.,;), : tangenual velocity component induced at the / panel control point by 
unit strength vorticity distribution on the wake panel of the first airfoil at time J,. 


e (B8'.,:), : normal velocity component induced at the / panel control point by unit 
strength vorticity distribution on the wake panel of the second airfoil at time 4,. 


e (B:.,;), tangential velocity component induced at the / panel control point by unit 
strength vorticity distribution on the wake panel of the second airfoil at time J,. 


@ = (4hp3,). : X Velocity component induced at the first airfoil wake panel midpoint 
with respect to the frozen local frame of reference by unit strength source distrib- 
ution on the /” panel at time /,. 


® (Akp;.), 2 ¥ velocity component induced at the first airfoil wake panel midpoint 
with respect to the frozen local frame of reference by unit strength source distrib- 
usion on the” panel at time 1,. 


@ (A%p,.), : X Velocity component induced at the second airfoil wake panel midpoint 
with respect to the frozen local frame of reference by unit strength source distrib- 
ution on the j” panel at time 4. 


6 It is important to note that the velocities at the midpoint of the wake panels do not include 
the effect due to the vorticity distributed along itself but it does include the contnbution due to the 
wake panel of the other airfoul. 


(A3-p3,), + ¥ Velocity component induced at the second airfoil wake panel midpoint 
with respect to the frozen local frame of reference by unit strength source distrib- 
ution on the j" panel at time £,. 


e =6(Bkp;,), + X velocity component induced at the first airfoil wake panel midpoint 
with respect to the frozen local frame of reference by unit strength vorticity dis- 
tribution on the j” panel at time /,. 


e =6(Bx>;,), + ¥ velocity component induced at the first airfoil wake panel midpoint 
with respect to the frozen local frame of reference by unit strength vorticity dis- 
tribution on the j” paneiamtumes. 


¢ (Bxpq;), 1X Velocity component induced at the second airfoil wake panel midpoint 
with respect to the frozen local frame of reference by unit strength vorticity dis- 
tribution on the j” panel at time 4,. 


© (Bx ps.), 1 ¥ Velocity component induced at the second airfoil wake panel midpoint 
with respect to the frozen local frame of reference by unit strength vorticity dis- 
tribution on the j* panel at time /,. 


e = (Aj,), : X velocity component induced at the A” core vortex with respect to the 
frozen local frame of reference by unit strength source distribution on the 
j* panel at times: 


e = (A}.), : ¥ velocity component induced at the h” core vortex with respect to the 
frozen local frame of reference by unit strength source distribution on the 
j* panel at time 4,. 


© 6(B;,), + xX velocity component induced at the h” core vortex with respect to the 
frozen local frame of reference by unit strength vorticity distribution on the 
j® panel at time J,. 


e =6(B.,), + ¥ velocity component induced at the h® core vortex with respect to the 
frozen local frame of reference by unit strength vorticity distribution on the 
jo manelwat Gime 


e (8 .,;), 1 x Velocity component induced at the vi meore vores o. Unit stremenm 
vorticity distribution on the wake panel of the first airfoil at time 4,. 


© (By .»;), 2 ¥ velocitv component induced at the A” core vortex by unit strength 
vorticity distribution on the wake panel of the first airfoil at time 4,. 


© (Bi yps), + X Velocity component induced at the h" core vortex by unit strength 
vorticity distribution on the wake panel of the second airfoil at time 4. 


(By yes), ° ¥ Velocity component induced at the A” core vortex by unit strength 
Vorticity distribution on the wake panel of the second airfoil at time 1,. 


2. New C’s Influence Coefficient 
The single airfoil definition is extended to include the effects of the second 


airfoil. 


e (C*,()), : normal velocity component induced at the 7” panel control point by 
unit strength +m” core vortex at time 4,. 
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(C;.(4)), : tangential velocity component induced at the / panel control point 
by unit strength om core vortex in the wake of the /* at time 1,. 


@ (Ck3,.(4), + X velocity component with respect to the frozen local frame of refer- 
ence induced at the first airfoil wake element by unit strength m" core vortex in the 
wake of the /* at time 4,. 


© (CrX»;,,(4), : ¥ velocity component with respect to the frozen local frame of refer- 
ence induced at the first airfoil wake element by unit strength m core vortex in the 
wake of the /” at time 4,. 


© (Crp (9), : X Velocity component with respect to the frozen local frame of refer- 
ence induced at the second airfoil wake element by unit strength mm core vortex in 
the wake of the /” at time 4,. 


® (Cr».-(9), : ¥ Velocity component with respect to the frozen local frame of refer- 
ence induced at the second airfoil wake element by unit strength m” core vortex in 
thie wweke of the 7” at timie {,. 


e (C7,()), : x velocity component with respect to the frozen local frame of reference 
mnauced at the h* core vortex by unit strength mm" core vortex in the wake of the /* 
at time t,. 


e (Ci,(J), : ¥ velocity component with respect to the frozen local frame of reference 
induced at the A” core vortex by unit strength mm” core vortex in the wake of the /” 
ac niet, 

3. Computation of the Time-Dependent Influence Coefficient 

As for the single airfoil , the influence coefficients of the wake element 
2223), Oa,.) avid (B:.,,), are computed as for the airfoil influence coefficients 
B: and B; with the subscripts NP3.NP4 replacing J. 

The x and y velocity components for the respective frozen local frames of ref- 
erence are obtained indirectly bv first calculating for the global X and Y velocity. These 
global velocities are then transformed to the frozen frames through a simple relationship 
bv the use of the respective angle of attack. It can be easilv shown that the velocity with 
respect to the frozen frame of reference is related to the global velocity by the following 


relationship. 


(Vink = [GV;,, cos a()], —[GV7,, sin a(D], 
(Vim = [GVim sin a()], +[GV;,, cos a()], (3.9) 


where V denotes a generic induced velocity with respect to the frozen frame and the 


precedent G denotes global velocities. 
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The global (GAfps,)a(GAxes Jal GA kira, )as «+: (GBhra,)(GA2).(GA!).(GB:), and 
(GB), are computed using equations 2.18 through 2.21] with 6, set to zero and subscript 


1 appropriately replaced. 
The global C’s coefficients will be computed identically as for the single airfoil 


case. The results are repeated here for the sake of completeness. 


Cosl(8;)e — (Orn 
(GCp(D)_ = — Ell (3.10) 
Sin{(@;), — (8,, 
(GChn(De = — Ei (3.11) 
where : 


(Finda = sf [(oert, = x2)? + Git, = Yn)? I 
Xi (eee 


vin, = 1/2 (¥, eed) 


X,t= xX Coordinate of mi core vortex at time #, 
y., = Y¥ coordinate of 71 core vortex at, times, 
“ plete 
, = tan Cay =a = ) 
hit, = Vip 
» aif Se es 
(6), = tan" xm x, Ve 


Also (GCkpsm())is (GCXpsm())as (GC Cipam Vas (GXpam Mae (GC), and (GC,()), are 
computed with equation 3.10 with 6, set to zero and the subscript 1 appropriately re- 


placed. 


E. NUMERICAL SOLUTION SCHEME 


1. The Flow Tangency Condition 
The unsteady flow tangency equation for the two airfoil system is a simple ex- 


tension of the single airfoil case. 


n 2n 


n 2n 
D> (Ag+ Y AGG del + 7Da> BY + VD). BF + (Paraame* Fl 
j=) j=n+) j=) j=n+!) 


a4 [yD] Bp vps), al [yw(2))Bivpa)x a5 
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R—1 
(Se Se(!))( at) ~ F(a = Dalia ONE (D —_ FO = Q 


m=] ni=1 
where 1= 1,2... 2n 
(Sele) 


where 1= 1,2... 2n and J%,,,,,, 18 evaluated by equation 3.1] at the # panel control point. 
Rearranging the equation we can express the source strengths explicitly as a 


function of the two vorticity strengths and a constant term as in equation 2.28 where : 


LHS. = > AN (gid + Satan (3.13) 


j=nt+) 
ims ed. S. = ay a AB; xp3)k Den (3.14) 
2n 
Second R.H.S. = hd Bay Bae > nd (3315) 
y j=n+) 
= = Sul i 
eeri.s. = —(() camettile — (1), -4( atk )( Bi wes) ~ 
SS(2) n 

He \ ~A(2)- \x( Bi xps)e 

k—| k-} 
q | DMCC nl . F(a - | LCi T ms _ FO (3.16) 

m=) m=] 


where 1 = 1,2... 2n 

For the code, the above is implemented in a 2N x (2+ 3) matrix in SUB- 
ROUTINE COEF and is subsequently solved by GAUSS elimination to obtain an ex- 
plicit expression for the source strength in terms of the vorticity distribution and a 


constant term as in the steady flow case as follows : 


Pat 


gq, = 61,¢ 9) > 02) <7 ee eel (3.17) 


where 


b1, : part of the source strength which depends on circulation y(1),. 


62, : part of the source strength which depends on circulation y(2),. 


63, : part of the source strength which is independent of the circulation. 


Equation 3.17 is critical in the simplification of the treatment for the 
non-linear coupled Kutta condition. 

2. The Kutta Condition 
For the unsteady flow, the Kutta condition can be written as the following : 


For the first airfoil 











(VoE- Wo =ALn- od ke = AS 
Rein (3.18) 
= 2SS(1) —— - = 
k k-} 
For the second airfoil 
, 2 rl a or pa) 
(lan), ~ (Va =a 2| a (Pop a Pn+) Ix = 2( - Vr 
yO) — (2 (3559) 
= 255(2) ————— 
Ep Tey 
The tangential velocity for the first panel can be written thus: 
2n n 2n 
y= - » (Brae + (1. 43 sea MO) yy Ai,+ 
j=!) c=) j=nt+i 
SSU) SS(2 
A(1) —— A; NPR R101) — y,(1)] za cit A; NPALY R—1(2) an y,(2)] a (Vv Say ef; tile + 
FY 'iciat DF nas = Fra(D) Ie + fecha DP ni(d FOI (3.20) 
m=) m=) 


Substituting equation 3.17 into 3.20, the following simplified relation for 1’; 1s obtained. 
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seein), +°2)7(2), + D3, (3.21) 


where 
2n n 
n 7. Sl 
a — » Bibl, +) 4y- — Ay.np3 
j=) j=) 
2n 2n 
n ny SZ 
> > Bie 5 DA iy Ay xps 
j=) j=nt!l 
: oa SSM 2» SS(2) 
D3, = - 2 Pb TAG) A1se3te-1) + ea) A(2) A; NPA k=] (Ze Gl V stream)! “ty ile + 
jJ= 
k=1 
+ px Chal F ny (l) - FO + AS 1CaF O- Oe (3.22) 
m=} Wi 


Sinular expressions can be obtained for the other trailing edge panels as follows : 


\ = D1,y()), + D2,y(2), + D3,, pus 
e mera (Mer 2 447 (2)y + D3, (3.24) 
J = = Di,,7(1), + D2>,74(2), + D3, (3.25) 


Where the coefficients D1, D2, D3 are obtained as per equation 3.22. We require the 


square of the tangential velocity for the Kutta condition 


(Vj)? =(D1,y(1), + D2,7(2), + D3,)’ 
= Dl?y¥(1)j + D2? (2), + D3? (3.26) 
+ 2D1,D2,7(1),7(2), + 2D1,D3,7(1), + 2D2,D3,(2), 


(Vi) =(D1,7(1), + D2,7(2), + D3,) 
Poel D2, (2), + D3. : (3.27) 
+ 2D}, D2,7(1)x7(2)p + 21, D3,y(1), + 2D2,D3,y(2), 


For the first airfoil equations 3.26 and 3.27 are then substituted into 3.18 to 
give, for the code: 


AAA(1)y(1) + BBB(1)y(2); + CCC()y(1), + 

DDD(1)y(2), + EEE(\)y(V)yy(2)e + FFF) == 0 (3.28) 
where 

AAA(1) = DN- DIP 


BBB(1) = D2; — D2; 


CCC(1) = 2D1, - D3, pl ape 
Ey ~ bpey 
DDD) = 2)D2.- 24 — D2. oem 
EEE) = 2(Dl,- D2, =o rere) 
SO Ep 
FFF(1) = D3? - D3 4 se > (3.29) 


nos 
A similar set of equation can be obtained for the second airfoil as follows: 

AAA(2) (1); + BBB(2)y(2)p + CCC(2)¥(1)y + 

DDD(2)y(2), + EEE(2)y(1)py(2), + FFFQ) = 0 EE!) 
Polen 

AAA(2) = DI?,,— Dl, 

BBB(2) = D22,, — D2, 

CCC(2) = 7 ON iy : DE — Dts, ¢ D32,) 


SS(2) 
DDD(2) = 2D2n4 * D3n41 — D2an * D3 an - Fp 


EEE(2) =) 21DY,, © D2, 9 lee 
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SS(2)4-)(2) 


PRR) = D3) un — D3 2 
Ey Eye y 


) (2.31) 
The above Kutta equations are non-linear and coupled. To linearize it, we adopt 
Newton's method, in which the initial iteration is linearized about the values of the 
previous time step: 


W(De = year + OE = DE = WCEP + 2v(Miclov(yy 
(De = y(pry + Sve = (CDEP = [y(Qy_p}? + 2v(2){716y(2)5 (3.32) 
where: 


SSE! sy yDR= rey 
b3(2YE<y(e) ss x22 = (Dey 


with x denoting the iteration counter. After dropping quadratic terms of dy and col- 


lecting like terms the Kutta condition for the first airfoil simplifies to the following: 
tema ))7(1),_, + CCC) + EEE(1)y(2),_, }oyQ)), + 

(250 8(1);(2),_, + DDD()) + EEE(I)y(1),_, }}oy(2), + 

- AAA(\)y(1),_, + BBB(1)7(2)2_, + CCC(1)y(1),_, + 


DDD(1);(2),_-, + EEE() (1) py (2)yp_-) + FFF(I) = 0 


(3-53) 
Similarly the linearised Kutta condition for the second airfoil can be written: 
{2AAA(2)y(1),_) + CCC(2) + EEE(2)y(2),_,}Ov(1)y + 
{2BBB(2)y(2),_, + DDD(2) + EEE(2)y(1),_-,}6r(2), + 
AAA(2)y(1);-, + BBB(2)y(2)p-) + CCC(2)W(Dg_-1 + 
DDD(2)y(2),_, + EEE(2)¥(1),_,y(2),_) + FFF(2) = 9 

(3.34) 


3} 


The above two linear equations with two unknowns ody(1), and 6y(2), can now be easilv 
solved with the same Gauss routine. The results are then back substituted into 
equations 3.32. This procedure is repeated until the corrections dy(1)s and 6+(2): are less 
than a prescribed tolerance. 
3. Convection of the Wake 
The resultant velocities of all core vortices are calculated using the frozen local 


frame of peter C This As resolved into components at the h” core vortex as follows: 
2n 


tine Plo (alk + ay (ae + 1p a 


j=nt+l j=n+] 
=F _ Vel Br ves + (yl Ba wpade 


k—] k-] 
tr jl halal 7 Fo a5 jl _ FO 
m=1 m=1 
mh m#h 
(3.35) 
2n 
(Vide = i (ale + >; [Ai(g)) iF ale v1), De. am y(2)x vm (Bray a (lp Fh 
j=n+) j=n+) 
a ae vp + (y(2)} (Ba, npa)ie 
+ {Yuet Ce AOA DN + {chon (Fn - FO 
m=) m=] 
mth mth 
(3.36) 


The location of the core vortices at the new time step is then computed with respect to 
the current time step frozen local frame of reference and then transformed to the new 
time step frozen local frame of reference. These coordinates are subsequently trans- 
formed to global coordinates so as to facilitate the calculation of the influence coeffi- 
cient. 
4. Disturbance Potential and Pressure Distribution 

The concept of disturbance potential introduced in the steady flow has plaved 
a Vital role in that it facilitates the synthesis of the flow field from simple flow field by 
simple superposition of the various singularities contributions. However, there has not 


been a requirement to solve directly for the disturbance potential. as our interest was 


on the disturbance induced velocity which is the spatial derivative of the disturbance 
potential. In our solutions, the concept of influence coefficient has allowed a direct 
evaluation of the disturbance velocity thus nullifving the requirements of obtaining the 
disturbance potential. 

The treatment for the unsteady flow, though it follows the same procedure as 
the steady flow case, viz the influence coefficient to obtain the disturbance velocity, still 
requires the disturbance potential for the computation of the pressure coefficient as can 


be seen in equation 3.4. Rewriting, we have: 


= sneer AE 2 (d;) (d;) = : 
(Cy) = Md = GE Me iH i= 1,2... 20 (3.37) 


where }.,... and }’,,,., are defined according to equations 3.1 and 3.3 respectively. 

From equation 3.37 we have written the rate of change of @ by a backward fi- 
nite difference approximation. This simplifies our iteration procedures tremendously as 
the @ from the previous time step exists at the current computation. The computation 
of the disturbance potenial @ is obtained through two steps. The difference in potential 
@ from upstream at infinity to the leading edge is computed and is then combined with 
the difference in the potential from the leading edge to the panel’s of interest control 
point. 

The present approach differs from the original single airfoil approach in that the 
velocity potential along the airfoil surface is computed via a finer grid using Gaussian 
quadrature’. This modification improves the results8 but the cost of doing it is a longer 
computation time. The definition of infinity has also been set to 100 chord lengths in 
comparison to 10 chord lengths used in the original code. However, the number of 
computation points and the first panel length from the leading edge upstream have not 
been changed. 

For completeness, the computation of the disturbance potential from the lead- 
ing edge to infinity is included here for both airfoils. It is essentially the same as the 
single airfoil case except that the disturbance potential is now computed relative to the 


global frame of reference instead of the airfoil fixed frame of reference as used in the 


7 Each panel is subdivided into 4 additional sub-panels and the tangential velocity is computed 
and integated over these srnaller panels with a weighting function to get the disturbance potential. 


8 One can check the improvement by computing the (5s) obtained through the difference be- 
tween the trailing edges and compares them with that obtained through the Kutta condition. 
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original code. We begin by selecting a straight line extending upstream in the direction 
parallel to V,.. The length of the line is set at 100 chord lengths. This line is divided into 
Z panels with the first element at the leading edge set equal to the single airfoil case for 
the purpose of ensuring that the panel size is comparable to the airfoil panel size. The 
panel size is then subsequently increased to take advantage of the inversely decaying 
induced velocities at larger distances. Using subscript f to denote these panel mid-points, 


we define the following influence coefficients: 


¢ (47,), : normal velocity component induced at the f* panel control point by unit 
strength source distribution on the J*panel at time 1,. 


(4}.), : tangential velocity component induced at the /” panel control point by unit 
stren oth source distribution on the j“panel at time 1,. 


e (87), : normal velocity component induced at the f” panel control point by unit 
strength VOrticity distribution on the J*panel at time 1,. 


e (8;,), : tangential velocity component at the /” panel control point by unit strength 
vorticity distribution on the "panel at time 4%. 


© =(Biyp3), + normal velocity component induced at the /” panel control point by unit 
strength vorticity distribution on the wake panel of the first airfoil at time 4,. 


¢ (Bixp;), + tangential velocity component induced at the f” panel control point by 
unit strength vorticity distribution on the wake panel of the first airfoil at time 4,. 


© (Brspa), : normal velocity component induced at the / panel control point by unit 
strength vorticity distribution on the wake panel of the second airfoil at time /,. 


e (B.-p,), tangential velocity component induced at the f” panel control point by unit 
strength vorticity distribution on the wake panel of the second airfoil at time 4. 


e (C7,()), +: normal velocity component induced at the / panel control point of the 
if“ DY unit streneth mm” core vores dummies. 


(C:,()), : tangential velocity component induced at the /" panel control point 
of the /* by unit strength mm” core wortex dt timer, 


The tangential velocity at the f* panel written as for the code is then: 
£ : p 


(Ve = - -) 50 e+ (1) Deir Scape 


j=nt+l 


Sse SS(2 
vr (Ap wp3)d?¢—1(1) — 7,1) + ane (Ay wea)id?x—r(2) — ¥(2)] + 


{Suc ‘im eT mal) — Fin(d) past {Hci Tipe — Trl i (3.38) 


m=1 
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valid for f=1,2 ...z. The disturbance potential at the airfoil leading edge is the sum of 
the products of the disturbance induced velocity at each panel and the panel length . 


(rede = - SK Vie (Apr — Ap] (3.39) 
ia 

The integral over the airfoil surface, as stated before, is now done over a finer 

grid. This requires the computation of the disturbance induced velocity over the smaller 

grids within each panels. Defining the total finer grid points in one panel by P9, we de- 


fine first the refined influence coefficient 
(Ai), 


tangential velocity induced at the # panel p” node due to unit strength 
source distribution on the j“ panel at time 2, 


The other refined influence coefficients have the same definition. The tangential com- 


ponent of the disturbance induced velocity at the #* panel m" node 1s then 


on n 2n 
(V4, an = > Bary (> Aly + y(2 Ve » Aip+ 
= j=) j=n+) 
SS(1 SS(2) 
+ A; “spa pl? pee) = gt) | + — A; vps pl? Ee a2) 
IY Chg) Final) — Fn(d) Ie + Yea ‘(OP ma() — Vn (0 fn (3.40) 
m=) m=) 


which is valid fori=1.2... 2n; p=1,2... P. Performing the line integration by summa- 


tion, the disturbance iin at the /* nodal point 1s then 


(Prode Ik = (Piedx + y de Vipetiiei Yo], forn> ie ite 


| (3.41) 
to] p 


= (Pedy >(%) ip)k kT} j+1 Wy) , for n>t2 le 


t=, p=) 


9 This includes the end points. 
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where y,,_, denotes the panel length and I”, the Gaussian quadrature weighting function. 





ter = i Mpa — A + Bpar — YP (3.42) 
Finally the disturbance potential at the /’ panel control point ts: 
(die - 1/2[ (Pnode dr Bi (Prode i+1)xl , i= i 1 On (3.43) 


5. Computation of Forces and Moments 
The C,, C, and C,, about the leading edge of each airfoil are calculated in exactly 
the same way as it is done for the steadv flow problem by integrating the pressure dis- 
tribution (See section D-4 of Chapter 2). 
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IV. DESCRIPTION OF COMPUTER CODE USPOTF2 


A. PROGRAM STRUCTURES AND CAPABILITIES 
1. Restrictions and Limitations 

The restrictions and limitations listed in the single airfoil documentation [Ref. 
1) still apply for the two airfoil case. Some efforts were made to optimise storage re- 
quirements versus computational repetitions; however as in the original code there is still 
much room for improvement. Krainer [Ref. 6 ] has improved the original code by 
combining subroutines and reducing repeated computation by introducing additional 
common variables. Some of his improvements have been added to this program. 

The computer system used is the Naval Postgraduate School IBM 3033AP. The 
current program fixes the maximum number of airfoil panels to 20010 and the maximum 
allowable time steps to 200. The computer currently has to be run with a minimum 
storage requirement of 2 Mbytes. A detailed computing time study for the program 1s 
not undertaken, for it not onlv changes with the number of nodal points selected but 
also with the time step increment as this has an effect on the rate of convergence for the 
wake panel iteration. An order of magnitude is given for one case run to give an ap- 
preciation of the computing time required. The system currently requires a total CPU 
run time of 200 seconds using an optimising compiler for a step input with a 0.025 
non-dimensionalised time increment for 26 time steps. 

At present. the program can run for two airfoils set at arbitrary distance and at 
different angles of attack undergoing anv of the following motions: 

1. In-phase and out-of-phase Step Input 

2. In-phase and out-of-phase Modified Ramp Input 

3. In-phase and out-of-phase Translational Harmonic Oscillation 
4. In-phase and out-of-phase Rotational Harmonic Oscillation 


5. Sharp Edge Gust Field Penetration 


2. Current Structures of USPOTF2 Main Program 
The flow logic for USPOTF2 (Unsteady Potential Flow for 2 airfoils) is illus- 


trated in Figure 6. It is essentially divided into three major modules, namely, the input 


10 These are combined total panels for both airfous. 
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- problem setup, Steady flow solution with its associated output and the Unsteady flow 


solution With its associated output. 


The first module includes subroutines INDATA, SETUP, BODY and NACA4S. 


The main functions of this module are: 


IL. 


Set up the problem formulation by reading in the necessary values and flag setting 
from filecode 1. 


If the flag fora NACA 4-digit or 5-digit of type 230XX is set, this module calls 
subroutines BODY and NACA45 to compute the local panel coordinates. 


If the flag fora NACA 4-digit or 5-digit of type 230XX is not set, this module then 
proceeds to read in the local panel coordinates from filecode 1. 


The local slopes and the airfoil perimeters are then calculated in preparation for the 
next module. 


The second module calculates the Steady flow solution. It calls subroutines 


NEWPOS, INFL, COEF, GAUSS, KUTTA, VELDIS and FANDM. The main func- 


tions of this module are: 


Ie 


Transform the local airfoil coordinates into global coordinates and compute the 
influence coefficients. 


Set up flow tangency equation as per equation 2.27 and solve for the source 
strengths (g,) as a function of the vorticity strengths [y(4] and a constant part where 
j}= 1,2... 2n and /= 1,2. 


Set up and solve the Kutta condition as per equations 2.25 and 2.26 for the 
vorticity strengths and back substitute to get the source strengths. 


Compute the total tangential velocity in the moving frame and compute the dis- 
turbance potenual (@,) and pressure coefficient [(C,)} QQ = 1,2... 2n). 


Finally, compute the aerodynamic coefficients of forces and moments -- 
CC 


This program can terminate in this module after all the steady flow parameters are 
obtained without necessarily running the Unsteady flow code. 


The third module calculates the Unsteady flow solution. It calls subroutines 


NEWPOS, INFL, COEF, GAUSS, KUTTA, TEWAK, PRESS, FANDM and 
CORVOR. The main functions of this module are: 


Ie 


For the particular unsteady motion, to compute the initial time step and new airfoil 
orientation with its associated local body velocities. 


Introduce the wake panel and assume an initial length and orientation to begin it- 
erations. 


Transform all local coordinates and panel slopes to global coordinates and global 
panel slopes. 
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4. Update influence coefficient and set up flow tangency equation. 


0) 


Solve for the source strengths in terms of the vorticity strengths and the constant 
part. 


6. Invoke the required Kutta condition and solve for the vorticity strengths!!; back 
Substitute to get the source strengths. 


7. Update the local velocities of the wake element and compute new length and ori- 
entation --- check for convergence. 


8. If wake element velocities have not converged, iterate with the new wake element 
geometries until convergencel2. 


9. Compute the total velocity, disturbance potential and the pressure coefficient. 
10. Compute the aerodynamic hft, drag and moment coefficient. 


11. Adjust time step either through an interative input!3 or bv a automatic time in- 
crement. 


12. Compute resultant local velocities of the core vortices. 


13. Convect the core vortices by the procedure described in Chapter III section E-3 


11 There are two possible solutions; only the vorticity distribution that ensures that the prod- 
uct of the relative tangential velocities of the upper and lower panels of the trailing edge is negative 
is accepted as solution. 


12 The convergence criterion is user specified through input data for TOL. 


13 This is accomplished by setting TADJ to be non-zero and is intended for use in conjunction 
with the viscous flow program. 
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READ INPUT FROM 
FILECODE 1 
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& PERIMETERS 


INPUT MODULE 


Figure 6. Flow Chart for USPOTF2 Computer Code. 


40 
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B. DESCRIPTION OF SUBROUTINES 
1. Subroutine BODY | 
This subroutine is called for the purpose of obtaining the local (x,y) coordinates 
of either a NACA XXXX or 230XX type airfoil. It is called by subroutine SETUP and 


it in turn calls subroutine NACA4S to obtain the airfoil thickness and camber distrib- 
utions. 
2. Subroutine COEF 
This subroutine was originally intended for the unsteady flow calculation for the 
single airfoil case. It is now modified to include the steady flow calculations. Its pur- 
pose is to set up the flow tangency matrix as in equation 2.27 for the steady flow and 
equations 3.13 through 3.16 for the unsteady flow case. These matrices are necessarily 
set up in this way so that the source strengths can be solved in terms of the vorticity 
Strengths and a constant by subroutine GAUSS as a linear system with three right hand 
sides. It is called bv the MAIN program. 
3. Subroutine COFISH (deleted for the two airfoil case) 
This subroutine was originally set up to serve the same function as subroutine 
COEF for the steadv case. It is now deleted for computational efficiency. 
4. Subroutine CORVOR 
This subroutine 1s called for the purpose of obtaining the convective velocities 
for all the wake core vortices with respect to the frozen local frame of reference of the 
current time step in accordance with equations 3.35 and 3.36 where all the influence 
coefficients are now locally computed (previously done in subroutine INFL) to save 
storage requirements since these are only required 1n this subroutine. The local influence 
coefficients are obtained indirectly bv first obtaining the global influence coefficients and 
transforming them viz equation 3.9. This subroutine is called by the main program 
nearing the end of the unsteady flow calculations before starting a new time step. 
5. Subroutine FANDM 
This subroutine is intended to calculate the overrall lift coefficient, drag coeffi- 
cient and moment coefficient about the leading edge for the two airfoils. In order to 
preserve the option of obtaining the x and y forces with respect to the respective local 
frames of reference, the original method rather than equations 2.34 through 2.36 1s im- 
plemented. This subroutine is called by the MAIN program in both the steady and un- 


steady flow computations immediately after computing the pressure coefficient. 
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6. Subroutine GAUSS 
This subroutine is the standard linear svstem solver that employs the well- 
known Gaussian elimination with partial pivoting and operates simultaneously on a user 
specified number of right-hand-sides. It is called by the MAIN program in both the 
steady and unsteady flow calculations. In order to use GAUSS, the coefficients of the 
augmented matrix must be set up so that GAUSS will return the solutions replacing the 
corresponding columns of the augmented matrix that were initially occupied by the 
right-hand-sides. The coefficient set-ups are done by subroutine COEF for both steady 
and unsteady flow problems. 
7. Subroutine INDATA 
This subroutine is intended to read in the first three to five cards of the input 
data depending on whether IFLAG = 0. The first three cards contain some description 
of the airfoil tvpe, problem definition, IFLAG information as well as the number of 
lower and upper panels. If IFLAG = 0, it will treat the airfoils as NACA tvpe airfoils 
and will proceed to read the NACA number and to calculate the thickness parameters 
that will be required by subroutine NACA45. This is the first subroutine called by the 
MAIN program. 
8. Subroutine INFL 
This subroutine generates most of the influence coefficients that are needed and 
shared by the different subroutines. It has been modified to include the steadv flow case 
for the purpose of reducing repeated computations. It utilises the known relative ge- 
ometrical parameters of the singularities to carry oul computation based on equations 
2.18 through 2.21 for the steady flow calculations and including 3.9 through 3.11 for the 
unsteady flow case. The MAIN program calls this subroutine in every iteration cycle 
of each time step so that the time dependent influence coefficients can be updated as 
and when neccessary. Time independent coefficients are computed once in the entire 
flow solutions]4.There are also time dependent influence coefficients that are independ- 
ent of the iterative cycle to obtain the wake panel orientation!5. Those influence coef- 
ficients involving the wake core vortices are also independent of the iteration cycle and 
are also updated once in each time step; but the process of the update is more compli- 
cated. This is due to the process that the wake need first to be convected with respect 


to the previous time step frozen frame of reference, transformed to the present frame of 


14 These are: on the airfoil panel by the panels on the same airfoul. 
15 These are: on the airfoil panel by the panels on the other airfoil. 
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reference and finally transformed to the global frame of reference by subroutine 
NEWPOS. Finally the influence coefficients involving the wake panels are calculated 
as frequently as the number of iterations take to find a converged solution. 
9. Subroutine KUTTA 
This subroutine is intended to solve the Kutta equation. It has been modified 
to include the steadv flow solution. After the source strengths have been determined by 
subroutines COEF and GAUSS in terms of the vorticity strengths and a constant, this 
subroutine invokes the Kutta condition as in equations 2.25, 2.26 for the steady flow 
case and 3.33, 3.34 for the unsteady flow case. The two linear equations!6 are then 
solved again with Gaussian elimination to obtain the vorticity distributions. For the 
unsteady case, we add in the additional requirement of finding the product of the 
tangential velocities at the upper and lower trailing edges. These are demanded to be 
negative as there are basically two possible solutions to the Kutta equations due to its 
original quadratic nature. 
10. Subroutine NACA45 
This subroutine is intended to calculate the camber and thickness distribution 
of the NACA 4-digit and the NACA 5-digit airfoils of type 230XX which share common 
thickness distributions with the 4-digit airfoil having the same thickness to chord ratio. 
This subroutine is called by subroutine BODY which is in turn called by subroutine 
SETUP and in turn called bv the MAIN program. 
11. Subroutine NEWPOS 
This is a new subroutine introduced as a result of setting up the five frames of 
reference. Most of the coordinates computation is done with respect to the respective 
local frame of reference. Its purpose then is to transform all coordinates in the respec- 
tive local frames of reference to the global frame of referencel7. This is necessary for the 
computation of the influence coefficients as well as in the convection of the wake core 
vortices. It facilitates the simple requirements of defining the airfoil once only with re- 
spect to its local frame of reference. Orientation and displacement of the airfoil in the 


two-dimensional plane is henceforth calculated through this subroutine. This subroutine 


16 For the unsteady case, the original two equations are non-linear and were subsequently 
linearised in the discussion in Chapter 3. 


17 It has a secondary function to obtain the slopes for the airfoil panels and the wake element 
panels. 
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is called from the MAIN program in several locations immediately after the local coor- 
dinates are computed. 
12. Subroutine PRESS 
This subroutine calculates the pressure distribution over the airfoil panels after 
the iterative solution for the unsteady flow problem has successfully met the convergence 
criterion. It first computes the tangential velocities at all panel control points using 
equation 3.40 (with p = — then performs the disturbance potential evaluation at 
the current time step according to equations 3.4] through 3.43. Together with the dis- 
turbance potential data obtained from the previous time step, it calculates the pressure 
distribution using equation 3.37. 
13. Subroutine SETUP 
This subroutine sets up the local panel nodal coordinates for MAIN program 
by reading the fourth through seventh data sets of the input file if IFLAG = 1 is set. 
It skips the data reading if IFLAG = 0 and proceeds to set up the node distribution and 
calls subroutine BODY to calculate the airfoil local coordinates. The node distribution 
adopts a cosine formula in order to have closely packed panels toward the leading and 
trailing edges for improvements in solution accuracy. Regardless of how the nodal co- 
ordinates are obtained , subroutine SETUP determines the local panel slopes and airfoil 
perimeter length. 
14. Subroutine TEWAK 
This subroutine is intended to calculate the resultant velocity components at the 
mid-point of the shed vorticity panel with respect to the current frozen frame of refer- 
ence. These velocity components are necessary to ensure that the correct shed vorticity 
panels’ length and orientation are established. This is the governing criterion for the it- 
erative solution scheme of the unsteady flow case. This subroutine is called by the 
MAIN program at every iteration cycle of each time step for the unsteady flow calcu- 
lation. 
15. Subroutine VELDIS 
This subroutine essentially performs the same function as subroutine PRESS for 
the steady flow case. It is thus redundant and could be deleted with some modification 
work necessary for subroutine PRESS. This work has been implemented by Krainer for 
the single airfoil case. This subroutine is called by the MAIN program in the steady flow 


computation. 
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C. INPUT DATA FOR PROGRAM USPOTF2 

The Input data is simular to that for the single airfoil case. Program USPOTF? reads 
its input data from filecode 1. An example of an input data file is attached in Appendix 
B for the case when the airfoil nodal coordinates are input by the user. Computer gen- 
erated airfoil coordinates are another option that can be selected if the airfoil chosen 
belongs to the family of NACA 4-digit or 5-digit airfoils of type 230XX. To do this, the 
IFLAG parameter is set to zero in the first item of the 4” set of data card and replace 
the 5" and 6” set of cards by two cards containing the particular airfoil NACA number 
for the two airfoils using format (15). Figure 7 contains an itemised description of the 


sequential input variables. 


D. OUTPUT DATA FROM PROGRAM USPOTF2 

The Output data is similar to that for the single airfoil case. Appendix C contains 
a sample output data generated by using the input data set from Appendix B. Due to 
the repetitive nature of the output as a function of time, onlv selective time set data are 
shown. The output data file begins with writing out what the program has read from 
the data file followed by the computed nodal coordinates only if thev are generated by 
the program: otherwise it proceeds to write the airfoil computed perimeter length. The 
next set of output data are the steadv flow solution parameters of distributed source 
strengths, vorticity strengths, pressure-velocity distribution, force-moment coefficient, 
potential at the control nodal points and the potential at the leading edge. In addition 
some output is given to allow an assessment of the accuracy of the flow solution. This 
consists of the normal component of the velocities at the panel control points and the 
integral of the disturbance tangential velocity along the airfoil contour. If only the 
steady flow solution 1s required, the output terminates here. 

For unsteady flow, in addition to the above, the unsteady flow parameters are also 
printed. This includes, at every time step, the iterative printout for the convergence of 
the length and orientation of the wake element panel, the unsteady solution parameters 
on the airfoil similar to that of the steady flow parameters as well as the trailing wake 
core vortices data. Again the same checking mechanisms are inserted on the airfoil 
panels. A detailed explanation on the output variable names 1s listed in Figure 8. All 


output data are non-dimensional quantities. 
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Data Set #1 
FD iLe 


Data Set #2 


TITLE 


Data Set #3 
NAIRFO 
XSHIFT 
Yonune 
Data Set #4 


IFLAG 


NLOWER 
NUPPER 


Data Set #5 
x@D)Gy(1) 


Data Set #6 
Viel) 


Figure 7. 


Format (I5) - 1 data card. 


- Number of title cards to be used in Data Set #2. 


Format (20A4) - ITITLE data cards. 
- Headings printed on output for case run 


identification. 


Format (15,2F10.6) - 1 data card. 

- Number of airfoil, in this case = 2. 

- Relative X distance of the 2 airfoil's pivot position 
with respect to the airfoil global coordinate system. 

- Relative Y distance of the 2 airfoil's pivot position 


with respect to the airfoil global coordinate system. 


Format (315) - 1 data card. 

- 0 if airfoil is NACA XXXX or 230XxX. 

- 1 otherwise. 

- Number of panels used on both airfoil lower surface. 


- Number of panels used on both airfoil upper surface. 


Format (6F10.6) if IFLAG = 1 - variable data cards. 

- local non-dimensionalised x-nodal followed by y-nodal 
coordinates for airfoil 1. Total of Nlowert+Nupper+l 
nodal points divided into 6 points per data card. 

Format (6F10.6) - variable data cards 

- local non-dimensionalised x-nodal followed by y-nodal 
coordinates for airfoil 2. Total of NlowertNuppertl 


nodal points divided into 6 points per data card. 


List of Input Variables. 
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Data Set #7 Format (6F10.6,13) - 1 data card 
ALPI(1) - Initial angle of attack (AOA) for airfoil 1 in degree. 
ALPI(2) - Initial angle of attack (AOA) for airfoil 2 in degree. 
DALP - Absolute change in AOA in degree for non oscillatory 
motions. 
Maximum amplitude of AOA change in degree for 
rotational harmomic motions. 
Non-dimensional rise time (V,t/c) of 
AOA for motion involving modified-ramp change in AOA. 
Non-dimensional oscillation (wc/V_.) 
for harmonic motions. 
The length from the leading edge to the pivot point 
for the local system. 
- Flag for in-phase and out-of-phase motion. 
= 0 oUt-of-phase m@tion. 
- 1 in-phase motion. 
Data Set #8 Format (8F10.6) - 1 data card. 
UGUST - Magnitude of non-dimensional gust velocity along 
global x-direction. 
VGUST - Magnitude of non-dimensional gust velocity along 


global y-direction. 


DELHX( 1) - Non-dimensional translational chordwise amplitude 


fOr airror, | 
DELHX( 2) - Non-dimensional translational chordwise amplitude 


for arretoil 2 





Figure 7 (cont’d) 
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Data Set #8 
DELHY( 1) 


DELHY( 2) 


PHASE( 1) 


PHASE( 2) 


Data Set #9 
TF 


DTS 


TOL 


TADJ 
SCLA 


NGIES 


Figure 7 (cont’d) 


(Cont'd) 


Non-dimensional translational transverse amplitude 
fOr alrtonie 

Non-dimensional translational transverse amplitude 
foOreairtoid 2 

Phase angle in degree between the chordwise and 
transverse translational oscillation with the latter 
as reference for the first airfoil. 

Phase angle in degree between the chordwise and 
transverse translational oscillation with the latter 


as reference for the second airfoil. 


Format (5F10.6,15) - 1 data card. 


Final non-dimensional time to terminate unsteady flow 
solution. 

Starting time step for non-osc. motions if TADJ = 0.0 
No. of computational steps per cycle for 

harmonic motion. 

Baseline time step size for all motions if TADJ ¥ 0.0 
Tolerance criterion for checking the convergence 
between successive iterations of (U,,), and Vike: 
Factor by which DTS will be adjusted. 

Steady lift coefficient for the single airfoil at the 
specified AOA 

Option for changing the unsteady Kutta condition to 
satisfy the tangential velocity as per Geising's case. 
QO equal pressure at the trailing edge panels. 


1 equal tgt velocities at the trailing edge panels. 
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TK 

TKM1 
ALPHA(L) 
OMEGA(L) 


UCL) 


VCL) 


NITR 

VXWC(L) 
VYWw(L) 
FAKE (CL) 


THETA( L) 


GAMK(L) 


X1(J) 
Wi) 
X(J) 
ml) 
QC J) 
ePCI) 
V(J) 


Figure 8. 


Bate sscep 27, 

Time step ¢;_ 7. 

Angle of attack of airfoil L at time ty. 

Rotational velocity (positive counter clockwise) at time 
GomtoraairtoilalL 

Chordwise translational velocity (positive forward) at 
tame ¢;, for airfoil L 

Transverse translational velocity (positive downward) at 
time t, for airfoil L 

Iteration number. 

lveragave solution of (U,), of airfoil L 

ieraGive Solucionsotale) of airfoil L 

Iterative solution of shed vorticity panel length A; 
Ofvairtoml L 

Iterative solution of shed vorticity panel orientation ©, 
Otgairrro1! I, 

Iterative solution of the strength of the current 
Vouticity distribution of airfoil L 

Panel number. 

local x-coordinate of the midpoint of jea panel. 

local y-coordinate of the midpoint of j£? panel. 

global X-coordinate of the midpoint of jen panel. 

global Y-coordinate of the midpoint of j®? panel. 
Strength of source distribution on the jen panel. 
Pressure coefficient at the midpoint of the jen panel. 
Total tangential velocity at the midpoint of the jen panel 


with respect to the moving local system. 


List of Output Variables. 
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VNC J) 


PHIK( J) 


BEC y ) 


INTGAMMA 


CD(L) 
CL(L) 
CMC L) 


M 
DUC) 


LEC) 


X1(M) 


LIC 


ZED) 


ec) 


X2(M) 


120} 


CURCCH iy) 


Normal velocity at the midpoint of the pea panel 


with respect 


to the moving local system. 


Potential at the mid-point of the jth panel at 


the current time step. 


Potential at the mid-point of the jth panel at 


previous time step. 


Integral of the disturbance velocity around the airfoil. 


Drag coefficient of airfoil L. 


Lift Coefficient of airiromuel, 


Pitching moment coefficient about leading edge 


Of aixrfiorl 


Trailing wake core vortex number. 


X-coordinate 
of airfoil 1 
Y-coordinate 
Of airfioi me! 
X-coordinate 
Of ~aizfoi) 71 
Y-coordinate 
Of airfoil! 
X-coordinate 
of airfiolie2 
Y-coordinate 
Of alrroiiaZ 
X-coordinate 
Of aim fore 2 
Y-coordinate 


Of eirftoil) 2 


of the center of the mf" core 
with respect to local system. 
of the center of the mf" core 
with respect to local system. 


of the center of the mf"? core 


with respect to global system. 


of the center of the mt” core 


with respect to global system. 


of the center of the m*? core 
with respect to local system. 
of the center of the m*" core 
with respect to local system. 
of the center of the mf" core 
with respect to global system. 
of the center of the mt" core 


with respect to global system. 


vortex 


vortex 


vortex 


vortex 


vortex 


vortex 


vortex 


vortex 


Circulation strength of the mt® core vortex of 


airfoil i 


Figure 8 (Cont’d) 
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V. RESULTS AND DISCUSSION OF CASE-RUNS 


USPOTF2 is primarily written as a follow-up to U2DIIF for the single airfoil. All 
the case-runs with the exception of the gust case will be presented. The approach in the 
gust case 1§ not consistent with the requirements for an irrotational flow field and will 
not be treated in this report. The step change in angle of attack (AOA) will be compared 
with Giesing’s for the same airfoil set at the final AOA, undergoing an impulsive start 
from rest. As there exist no comparison data for the other sub-cases, the results will thus 
not be as extensive as the step input. They are documented here for the sole purpose 


of illustrating the capability of the code. 


A. STEP CHANGE IN ANGLE OF ATTACK 
I. Case Run Definition 
Consider a first case of two airfoils initially at zero AOA to the free stream J’, 
which undergo an out-of-phase step change in AOA (a,,,,) at time %. Consider a second 


case of two airfoils at rest, set intually at an AOA of a,,,, out-of-phase with one another, 


step 


given an impulsive start to I’_. The above two cases are equivalent within the thin 
airfoil approximation. Case 1 is computed bv USPOTF2 and case 2 1s obtained by 
Giesing’s computational analvsis. 
2. Differences between USPOTF2 and Giesing’s code 
While both codes use an extension of the PANEL method to solve for the un- 
steady potential flow solution for two airfoils, perfect correlation is not possible for se- 
veral reasons. The reader is referred to References [5] and [7] for a detailed description 


of Giesing’s approach. Some basic differences are: 


e USPOTF2 models the wake vortex sheet by a wake element and a series of point 
vortices shed through the wake element. This follows the approach of Basu and 
Hancock with the necessary assumptions on the wake element characteristics. 
Giesing treated the wake vortex sheet to comprise of a distributed line vortex 
which is convected at the local fluid velocity at the vortex location assuming that 
the small portion of the wake being shed does not contribute to the convection of 
the wake. 


e The time step increment in USPOTF2 is a parameter that affects the overrall results 
of the unsteady flow in that by having small time step. the point vortices being shed 
will be greater in number but weaker in strength, while for bigger time step, the 
point vortices become fewer but have effectively stronger strengths. The time step 
increment in Giesing’s code is important in that the assumption that the small 


Be 


portion of the wake being shed does not contribute to the convection of the wake 
is only exactly true in the linut when the time step becomes Zero. 


Giesing uses the Adam's formula of varving degreel8 to convect the wake vortex 
sheet while USPOTF2 uses only a predictor algorithm. 


The treatments of the circulation I'(/J) are different for both codes. USPOTF2 ex- 
tended the influence coefficient concept to the unsteady flow regime and solved for 
the circulation using the Kutta condition of equations 3.33 and 3.34. All calcu- 
lations Were done in the moving frame in the presence of the unsteady wake for- 
mations. Giesing considered the circulation to be a combination of the 
quasi-steady circulation and the circulation due to the vortex wakes. From Refer- 
Sivas || 3) 


Pr) = PoA0+ TY) Gan 


Where I,(/ 1s the circulation required to satisfy the Kutta condition on body (/) as 
it moves through the fluid when it 1s assumed that the bodv does not shed anv 
vorticity and T(/) 1s that circulation required to satisfv the Kutta condition on body 
(/) as it travels through the flow field generated by the vortex wakes shed bv the two 
bodies and the flow field generated by the circulatory flow about the other bodv. 


Giesing’s Kutta condition reqires equal tangential velocities at the upper and lower 
surface panels at the trailing edges while the Kutta condition of USPOTF2 pre- 
scribes equal pressures. In order to have a meaningful comparison, USPOTF2 in- 
cludes an option for equal tangential velocities. 


USPOTF2 requires the actual computation of the total velocity potential from a 
combination of an onset flow potential, disturbance potential and an assumption 
of a reference potential. Giesing uses the technique of the Douglas Potential Flow 
program Which treats the potential as the combination of the quasi-steady potential 
and the potential due to the Vortex wakes where these two terms are defined as 
before. The total velocity potential. in his case. need not be computed but only the 
time derivative of it Which is then written 1n terms of previously Known parameters. 


The number of panels and its distribution on the airfoils are different for both 
codes. This will be accentuated at the trailing edges and will lend itself to differ- 
ences in the trailing edge panel distribution. 

3. Results and Discussions 


Figure 9 shows Giesing’s results for the Von Mises 8.4 per-cent thick, symmet- 


rical airfoil undergoing an impulsive start. This can be compared with Figure 10 which 


shows the result obtained with USPOTF2 for the same time step with NGIES set to 


onel9. The results are not perfectly correlated but the quality and order of magnitude 


agreement are excellent. Figure 11] gives essentially the same result but with NGIES set 


18 This is essentially a predictor-corrector algonthm. 


19 This Kutta condition results in equal tangential velocities at the trailing edge panels of the 


airfouls. 
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to zero20. Surprisingly, the results obtained for both tvpes of Kutta condition turn out 
to be quite simular. This is seen especially in the pressure coefficient plots when the two 
plots are put together as seen in Figure 12. The pressure coefficient agrees over 70 per- 
cent of chord length at the lower surface and over 90 per-cent of chord length at the 
upper surface with the greatest discrepancy at the trailing edge. Further comparison for 
smaller angles of attack would be of interest to see whether this similarity is generally 
true. 

Figure 13 compares the aerodynamic characteristics when the vertical distance 
YSHIFT is varied. Figure 14 gives the time variation for a larger time step of the 


normalised lift, moment coefficient and drag coefficient for the particular case of 
YSHIFT = 2.0. 


20 This Kutta condition results in equal pressure coefficients at the trailing edge panels of the 
airfouls. 
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(a) Pressure Distribution and Vortex locations. 


Figure 9. Giesing’s Calculated Results: Impulsive Start for an 8.4 per-cent thick 
Von Mises airfoil set ata = 0.8 radians for YSHIFT = 2.0 [reproduced 


with permission from Ref. 5]. 
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_ (b) Time History of the Lift Coefficients. 


Fieure 9 (Cont'd) 
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(a) Pressure Distribution at t} "o/c 


Figure 10. 





USPOTF2 Results obtained with USPOTF2 for equal velocities at the 

trailing edge.: Step change in AOA for an &.4 per-cent thick Von 
Mises airfoils placed at 2 chord length vertical distance with iniual 
AOA = 0.0 radians and final AOA = 0.8 radians pivoting at the 
leading edges. 


5§ 


(b) Wake pattern at tVoo/c = 0.65 





Figure 10 (Cont’d) 
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0.2 0.4 0.6 
NON—DIMENSIONALISED TIME 


(c) Time History of the Lift Coefficients. 





Figure 10 (Cont’d) 
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(a) Pressure Distribution at t}_/c 


Figure 11. 





USPOTF2 Results obtained with USPOTF2 for equal pressures at the 
trailing edge.: Step change in AOA for an 8.4 per-cent thick Von 
Mises airfoils placed at 2 chord length vertical distance with initial 
AOA = 0.0 radians and final AOA = 0.8 radians pivoting at the 


leading edges. 


61 


(b) Wake pattern at tPoo/c = 0.65 





Figure 11 (Cont’d) 


0.2 0.4 0.6 
NON—DIMENSIONALISED TIME 


(c) Time History of the Lift Coefficients. 





Figure 11 (Cont’d) 
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NGIES = 0: Equal Pressure Coefficient at trailing edge panels. 


NGIES = 1: Equal Tangential Velocities at trailing edge panels 





Figure 12. Pressure Distributions with different Kutta condition. 


i——e— SINGLE AIRFOIL 
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i—-e-- YSHIFT = 3 
ooo) YSHIFT a 5 
poote- YSHIFT = :25 


(a) Pressure Distribution at r}/c = 0.65 





Figure 13. Pressure Distribution and Lift History as a function of the Vertical dis- 


tance between the airfoils. 
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(b) Time History of the Lift Coefficients. 


Figure 13 (Cont’d) 
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(a) Time History of the Lift Coefficients. 





Figure 14. Step Change in Angle Of Attack: Step change in AOA for an 8.4 
per-cent thick Von Mises airfoils placed at 2 chord length vertical] dis- 
tance with initial AOA = 0.0 radians and final AOA = 0.8 radians 


pivoting at the leading edges. 
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NON-DIMENSIONALISED TIME 


(b) Time History of the Moment Coefficients. 





Figure 14 (Cont’d) 
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(c) Time History of the Drag Coefficients. 





Figure 14 (Cont’d) 
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(d) Time History of the Circulation 





Figure 14 (Cont’d) 


B. OTHER SUB-CASES 

The reader is referred to the results and disussions of Reference [ 1 ] for the com- 
parison of the single airfoil case with existing codes in the phase and magnitude re- 
lationship of the aerodynamic coefficients to the forcing function as well as the wake 
convection. The present results show that the same trend is maintained as in the single 
airfoil case with no significant phase shift but with a general magnitude change due to 


an effective ground effect caused by the presence of the second airfoil. 
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1. Mfodified Ramp Change in Angle of Attack 


As for the single airfoil, the modified ramp is defined mathematically as follows: 


0 La) 
a(t) = d6a(3—2/t)7/e? O< 1 <t (5.2) 
Oo {> Tt 


Where da is the magnitude of the AOA change and 7 is the rise time for the AOA to 
reach its final value. Figure 15 treats the case of the 2 airfoils undergoing an out-of- 
phase modified ramp change in the angle of attack of 0.1 radians with a time constant 
of 1.5. The plots show the effects on the aerodynamic coefficients due to the presence 
of the second airfoil. The corresponding aerodynamic coefficients of the single airfoil 
are plotted for comparison purposes. 
2. Translational Harmonic Notion 

The code is capable of computing the unsteady flow solution for any general 
translational motion described by a chordwise and a transverse component bearing a 
given phase relationship with the restriction that the 2 airfoils move onlv in-phase or 


out-of-phase. The translational harmonic motion is described by 


h, (2) 


: 0h, Sin (a2) 


h(t) = dh, Sin (wt +4) (5.3) 


where qw 1s the oscillation frequency, / 1s the phase angle between the chordwise and 
transverse oscillation and of, and of, are the magnitudes of chordwise and transverse 
oscillations respectively. The case-run considered in this section relates to a pure heav- 
ing or plunging motion. A NACA-0015 airfoil is chosen for the case-run. The airfoils 
are set at zero radian angles of attack and subsequently given an out-of- phase plunging 
oscillation at an amplitude of 0.018 chord length at a non-dimensionalised frequency of 
1. Figure 16 shows the aerodynamic coefficients for the two airfoils and compares them 
with the single airfoil case where applicable. Note that the plots for the trailing wake 
uses different scales for the xX andy axes. 
3. Rotational Harmonic Motion 
The treatment of the harmonic pitching motion 1s similar to the modified ramp 


case. As for the other sub-cases, the airfoils are restricted to in-phase and out-of phase 


yy 


motion. The case of the out-of-phase motion will be treated here. The harmonic 


pitching oscillation 1s described by: 
a(t) = da Sin (ws) (5.4) 


where 6a and w are the amplitude and frequency of the harmonic oscillation respec- 
tively. Figure 17 shows the results of the 8.4% thick Von Mises symmetric airfoil os- 
cillating at an amplitude of 0.1 radian at a reduced frequency of wc/V,, = 20.0 about the 
leading edge. Again the plots are given together with the single airfoil case undergoing 


the same motion for comparison purposes. 
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(a) Time History of the Lift Coefficients. 





Figure 15. Niodified Ramp Change in Angle Of Attack: Ramp change in AOA 
for an 8.4 per-cent thick Von Mises airfoils placed at 1 chord length 
vertical distance with initial AOA = 0.0 radians and final AOA = 0.1 


radians, rise time of 1.5 pivoting at the mid chord. 
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(b) Time History of the Moment Coefficients. 





Figure 15 (Cont’d) 
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(c) Time History of the Drag Coefficients. 





Figure 15 (Cont’d) 
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(d) Time History of the Circulation 


Figure 15 (Cont’d) 
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(e) Pressure Distribution at rt} o/c = 2.3 





Figure 15 (Cont’d) 
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Figure 15 (Cont’d) 


78 


NON—DIMENSIONALISED TIME 


(a) Time Historv of the Lift Coefficients. 





Figure 16. Harmonic Plunging Motion: Translational harmonic plunging AOA 
for NACA-0015 airfoils placed at 1 chord length vertical distance set 
at AOA = 0.0 radians with plunging amplitude of 0.018 chord length 


at a reduced frequency of 1. 
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(b) Time History of the Moment Coefficients. 





Figure 16 (Cont’d) 
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(c) Time History of the Drag Coefficients. 





Figure 16 (Cont’d) 
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(d) Pressure Distribution at tr} j/c = 12.6 





Figure 16 (Cont’d) 
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Note : Different scales for x and v axes 
(e) Wake Pattern at rt} /c = 12.6 





Figure 16 (Cont’d) 
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NON—DIMENSIONALISED TIME 


(a) Time History of the Lift Coefficients. 





Figure 17. Harmonic Pitching Motion: Rotational harmonic pitching AOA for 
NACA-0015 airfoils placed at 1 chord length vertical distance set at 
AOA = 0.0 radians with pitching amplitude of 0.1] radian at a reduced 


frequency of 4 pivoting about the leading edges. 
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(b) Time History of the Moment Coefficients. 





Figure 17 (Cont’d) 


85 


NON—DIMENSIONALISED TIME 


(c) time History of the Drag Coefficients. 





Figure 17 (Cont’d) 
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1 
NON—DIMENSIONALISED TIME 


(d) Time History of the Circulation. 





Figure 17 (Cont’d) 
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(e) Pressure Coefficient at tr, 





Figure 17 (Cont’d) 
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Figure 17 (Cont’d) 
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VI. CONCLUSION 


A. GENERAL COMMENTS 

USPOTF2 has been developed as an intermediate stage to getting a solution for a 
full cascade undergoing unsteady motion. In itself, it can be used to simulate unsteady 
wing--tail, wing--canard, wing--flap, wing--aileron, horizontal tail--elevator, vertical 
tail--rudder and a whole host of moving--stationary airfoil interactions. These simu- 
lations require a little addition to the main program and would appear as modules 
modelling the variation of the global angles of attack and relative displacements as per 
all the sub-cases done. The subroutines would not be affected by the severity of the test 
case except for possibly subroutine NEWPOS. 

Validation of USPOTF2 has been done against Giesing’s code for the particular case 
of a step change in angle of attack. This was shown to have good correlation. However, 
it can not be said that the agreement will hold for other angles of attack and more sub- 
case runs for smaller angles of attack should be made. In the same manner, the success 
of the step change in AOA in no way validates the other sub-cases which, for the time 


being. remain unproven. 


B. ENHANCING USPOTF2 PROGRAM CAPABILITY 

As noted above, USPOTF2 needs to be tried more extensively either with the exist- 
ing sub-cases or with new sub-cases against existing numerical or experimental results. 

The software in itself has an implicit weakness in the numerical computation of the 
velocity potential. The velocity potential at the leading edge is obtained by integrating 
the velocity field from the leading edge to a point 100 chord lengths upstream of the 
leading edge. Two assumptions are made: First, the disturbance velocity at the point 
100 chord lengths of the leading edge approaches zero. Second, the contribution to the 
velocity potential due to the integration from the point 100 chord lengths upstream of 
the leading edge to upstream infinity must not change in time. The coefficients of pres- 
sure Will be accurate only, if those two assumptions hold (at least in an approximate 
sense). This however does not give a very ‘exact’ solution as by increasing the chord 
length by 900 per-cent to 1000 chord lengths, there is a change in the pressure coefficient 
at the trailing edge of about 8 per-cent for the first time step. Krainer has implemented 
an analytical solution to the potential problem for the single airfoil. A similar procedure 
to upgrade USPOTF2 should be considered. 
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Other improvements to the code would be the reduction of redundant computations. 
An obvious example would be the subroutines VELDIS and PRESS which essenually 
do the same work for the steady and unsteady flow respectively. Again Krainer [Ref. 
6] has improved the original code for the single airfoil (U2DIIF) and though some of 
his improvements were implemented in USPOTF2, there still remains a task to do the 
full job completely. 

Finally, the original primary objective of extending the code to solve for the un- 
steady flow solution for 3,4 airfoils and leading to a full cascade still remains a big task, 
not just for the programmer but also for the computer system in terms of computational 


storage and time requirements. 
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APPENDIX A. USPOTF2 SOURCE LISTINGS 
ccccccccccccccccccccccccccccccccceccecccccccceccccccccecccccccecccccccccce 


C C 
E PROGRAM USPOTF2 E 
C VERSION 1 c 
e SEPTEMBER 88 a 
C é 
c UNSTEADY MOTION FOR TWO AIRFOILS IN POTENTIAL C 
c INCOMPRESSIBLE FLOW C 
E USING PANEL METHODS BASED ON THE HESS & SMITH € 
@ AND WAKE MODEL BASED ON BASU & HANCOCK C 
C C 
ccccccceccccecececeececcceececececceccecccceceeccceuceccescecc rece eeeaan 
COMMON /BOD/ IFLAG,NLOWER,NUPPER ,NODTOT, X( 202) ,Y( 202), 

fs COSTHE( 201) ,SINTHE( 201) ,SS(2),NP1,NP2,NP3,NP4, 

rs NPS ,XSHIFT, YSHIFT ,NAIRFO, X1(202) ,YI(202), 

# COSTHL( 201) , SINTHL( 201) 


COMMON /WAK/ VYW(2),VXW(2),WAKE(2) ,DT 

COMMON /WAK2/ VYWK(2),VXWK(2) 

COMMON /SING/ Q(200) ,GAMMA( 2) ,QK(200) ,GAMK(2) 

COMMON /CORV/ CV(2,200),XC(2,200),YC(2,200) ,M,TD,CVVX(2, 200), 
+ CVVY(2,200),XCI( 2,200) ,YCI(2,200) 

COMMON /POT/ PHI( 200) , PHIK( 200) 

COMMON /GUST/ UG( 200) ,VG( 200) ,XGF,UGUST, VGUST 

COMMON /EXTV/ UE( 200) ,VN(200) 

COMMON /PARD/ NACAD(2),TAUD(2) ,EPSMAD(2) , PTMAXD(2) 

COMMON /COF/ A(201,211),KEQNS 

COMMON /GEOM/ SINALF(2),COSALF(2),OMEGA(2),UX(2),UY(2), PIVOT, 
+ XPRM , YPRM 

COMMON /CPD/ CP(200),SCL,T,SCM,SGAM 

COMMON /GIES/ NGIES 

DIMENSION XXC(2,200),YYC(2,200),TOL1(2) , TOL2(2),THENP1(2), 
+SINANG( 2) ,COSANG(2) , COSDA(2) , SINDA(2) , DHX(2) , DHY(2) , ALP(2), 
+ALPHA( 2) ,ALPI( 2) ,ANGLE(2) , DELHX(2) , DELHY( 2) , PHASE( 2), PHA(2), 
+HX(2) , HXO( 2) ,HY(2) ,HYO(2) 


C 
C INITIALISATION 
C 
Jel = 3. 1415926585 
M = 0 
XPRM = 0.0 
day = 0.0 
C 
C INPUT FROM FILE CODE 5 AND SET UP LOCAL PANEL NODES AND SLOPES 
C 


WRITE (6,1003) 
1003 FORMAT (/////,' DATA READ FROM FILE CODE 1',//) 
CALL INDATA 
CALL SETUP 
READ (1,502) ALPI(1),ALPI(2) ,DALP,TCON, FREQ, PIVOT, IPHASE 
WRITE (6,502) ALPI(1),ALPI(2),DALP,TCON, FREQ, PIVOT, IPHASE 
501 FORMAT (8F10. 6) 


qNdaqmgc. 


502 FORMAT (6F10. 6,13) 
READ (1,501) UGUST,VGUST,DELHX( 1) ,DELHX(2),DELHY(1),DELHY(2), 


+ PHASE( 1), PHASE( 2) 
WRITE (6,501) UGUST, VGUST, DELHX( 1) , DELHX( 2) , DELHY( 1) , DELHY( 2), 
+ PHASEC)) PHASE(2) 


READ (1,503) TF,DTS,TOL,TADJ,SCL,SCM,SGAM,NGIES 
WRITE (6,503) TF,DTS,TOL,TADJ,SCL,SCM,SGAM,NGIES 
SQ FORMAT (7F10.6,13) 
IF (IFLAG .EQ. 0) WRITE (6,1005) 
1005 FORMAT (///,' COORDINATES OF AIRFOIL NODES’, 
+ //,3X, X/C’ ,6X,' Y/C’,/) 
Meee 60. 0) White (6,1010) (XIC 1}, wit 1), 1=i 2*CNDDTOT+1)) 
1010 FORMAT (F10.6,F10. 6) 
WRITE (6,1000) NAIRFO 
1000 FORMAT(//,' TOTAL NO OF AIRFOILS = ',14,/) 


STEADY FLOW CALCULATION AT ALPI(I) 
INITIAL DEFINITION 

NP 1 = NODTOT + 1 

NP2 =52,> NP] +l 

NP3 = NAIRFO*NODTOT#+1 

NPS = NP3+1 

NP5 = NP4+1 


DO 101 I = 1,2 
ALPHA(I) = ALPICI) 
IF (ALPHA(I) .GT. 90.) GO TO 200 
COSALF(I) = COS(ALPHA(I)*PI/180. ) 
VOU SMOALF(I) = SIN(ALPHA(1I)*PI/180. ) 
CALL NEWPOS(0) 
DO 1100 L = 1,NAIRFO 
1700 WRITE (6,1020) L,SS(L) 
WRITE (6,1040) XSHIFT,YSHIFT 
1020 FORMAT(//,°§ AIRFOIL(' ,12,°) PERIMETER LENGTH = ',F10.6,/) 
Meoo Former (//,' XSHIFT = ',FS.1,' YSHIFT = ',F5.1,/) 
Dommo2 I = 1,2 
102 WRITE (6,1030) I,ALPHA(1) 
1030 FORMAT (//,' STEADY FLOW SOLUTION AT ALPHA(',I2,') = ',F10.6,/) 
IF (SCL. NE.0) WRITE (8,1035) 
dese bom 703%,’ DIME”,4%,'GL(1)/SCL' ,1X,'CL(2)/SCL' ,1X,*CM(1)/SCM', 
ee 602) /SCM ,3X, CD(1)' ,4X,'CDC2)' ,4X, 
+'GAMK(1)/SG' ,1X,'GAMK(2)/SG' ,//) 
[ee SeimEQ.0) WRITE (8, 1036) 
1036 HomMAgcsx ,'TINE' ,6X,'CL(1)',5xX,'CL(2)' ,5X, CM(1)°, 
Foe (2) ,SX, CD(1)° ,4X, 'CD(2)'//) 
CALL INFL (0) 
GALL COEF (0) 
CALL GAUSS(3,0,0) 
CALL KUTTA(NITR, PVTAG) 


CALL VELDIS 

SIN1 = SINALF(1) 
SIN2 = SINALF(2) 
COS1 = COSALF(1) 
COS2 = COSALF(2) 


a5 


De & | 


SB iat i 


100 


CALL FANDM(SIN1,SIN2,COS1,COS2) 


INITIALISATION FOR UNSTEADY FLOW CALCULATION TO BEGIN 


DA = 0.0 

XGF = 0.0 

DO 100 L = 1,NAIRFO 

HX(L) = 0.0 

HY(L) = 0.0 

HXO(L) = 0.0 

HYO(L) = 0.0 

ANGLE(L) = ALPI(L)*PI/180. + ATAN(VGUST/(1.+UGUST) ) 
OMEGA(L) = 0.0 

ALP(L) = ALPI(L) 
DHX(L) = 0.0 

DHY(L) = 0.0 

COSANG(L) = COS(ANGLE(L)) 
SINANG(L) = SIN(ANGLE(L)) 
UX(L) = 0.0 

UY(L) = 0.0 

COSDA(L) = 1.0 

SINDA(L) = 0.0 

KIG = (L -1)*NODTOT 
VXW(L) = COSALF(L) 

VYW(L) = SINALF(L) 

GAMK( L) = GAMMA(L) 


PHAC(L) = PHASE(L)*PI/180. 
DO 100° 1G =i. hepion 


UGCIGTKIG ae — O20 
VGC IG+KIG) = 0.0 
A) = 0.0 
TOLD = 0.0 


RIGID BODY MOTIONS OF AIRFOIL 


IF (FREQ .NE. 0.0) GO TO 
IF (DALP .EQ. 0.0) GO TO 
IF (TCON .NE. 0.0) GO TO 
IF (IPHASE .EQ. 0) GO TO 
ALPHA(1)= ALPI(1) + DALP 
ALPHA(2)= ALPI(2) + DALP 
GO TO 5 

ALPHA(1)= ALPI(1) + DALP 
ALPHA(2)= ALPI(2) - DALP 
po 6 I = 1,2 

COSALF(I) = COS(ALPHA(1)*PI/180. ) 

SINALF(I) = SIN(ALPHA(1)*PI/180. ) 

CONTINUE 

CALL NEWPOS(0) 

DT = DTS 

TD = DTS 

GO TO 60 

IF ((UGUST .EQ. 0.0) . AND. (VGUST .E0. 0.0)) GO TO 200 
DT = DTS 

TD = DTS 

GO TO 60 


fF WH rH 
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qaqa 


oo qag 


aan 


i DL wa2. OPE) (FREQ. DIS) 
TD = DT 
60 =f = DT 


Mrerre (6,1051) 
1051 FORMAT Gee ee ee Fereaie se ded dese ae ae ak ae ae ae ae ae ae ae ve se ae dese ae ae ae ee eM eae A eI HII | 


+ ' %%* BEGIN UNSTEADY FLOW SOLUTION Wied! 7 
a * Heve dese se dese ve ve vere Heese ae He Hes He ae TE NE Ie Te Te eI Me IH IIH HA TEI | + 


URRY = 0 
40 M =M+1 
im ({T .GT. TE) GO TO 200 
STORE CORE VORTEX COORDINATES FOR TIME STEP ADJUSTMENTS 
fgech . BO. 1) GO TO 50 
er 51 L = 1,NAIRFO 
Der 51 I = 1,M-1 
meeo(L,1) = XCI(L,I) 
DS) mro(L,1) «= YCICL,I) 


50 i (FREQ .NE. 0.0) GO TO 11 
me (OAEP .E@. 0.0) GO TO 22 
PeeCrCGN .NE. 0.0) GO TO 33 


STEP CHANGE IN AOA 
TE Smaps — NE. 0.0) GO TO 70 
TE MReREMENTAL PROGRESSIVE TIME STEP IS REQUIRED 
Wee. 1D — FLOATCMF1)*DTS ... OTHERWISE CONSTANT TIME STEP 
Tp = DTS 
GC, 19570 
MODIFIED RAMP CHANGE IN AOA 
5 me (1 .GT. TCON) GO TO 34 
DAL =o eee oa 2. /TCON) *( 1/ ICON) 222 
OMEGA( 1) = - (DALP*PI/180.) * (6.*T/(TCON*TCON)) * (1. -T/TCON) 


me CIPHASE .EQ. 0) GO TO 41 
ALPHA(1)= ALPI(1) + DAL 
ALPHA(2)= ALPI(2) + DAL 
OMEGA( 2) = OMEGA(1) 
GO TO 42 

41 ALPHA(1)= ALPI(1) + DAL 
ALPHA(2)= ALPI(2) - DAL 


OMEGA( 2) = - OMEGA(1) 
G2, bores Tia 1,2 
COSALF(I) COS(CALPHA(I)*PI/180. ) 


SINALF(T) = SINC ALPHA(CI)*PI/180. ) 


DA = ALPHA(I) - ALP(I) 
COSDA(I) = COS(DA*PI/180. ) 
SINDACT) = SINCDA*PI/180. ) 
DHX( I) = PIVOT * (1. -COSDA(I)) 
DHY( I) = - PIVOT * SINDA(TI) 
UY( 1) = PIVOT * OMEGA(1) 

43 CONTINUE 
MICON =4M 
CALL NEWPOS(0) 
GO TO 70 


aQqQaannqn 


34 


45 
46 


44 


Haka 


120 


Za 


LO 


DAL = 0.0 

IF (CIPHASE 2 EO Opec e 1eeas 
ALPHA(1)= ALPI(1) + DALP 
ALPHA(2)= ALPI(2) + DALP 


GO TO 46 

ALPHA(1) = ALPI(1) + DALP 
ALPHA(2) = ALPI(2) - DALP 

DO 44 I = 1,2 

COSALF(I) = COS(ALPHA(1)*PI/180. ) 
SINALF(I) = SIN(ALPHA(I)*PI/180. ) 
DA = 0.0 

COSDA(I) = 1.0 

SINDA(I) = 0.0 

OMEGA(I) = 0.0 

DHX(I) = 0.0 

DHY(1) = 0.0 

UX(1) = 0.0 

UY(I) = 0.0 

CONTINUE 


CALL NEWPOS(0) 

IF (TADJ .NE. 0.0) GO TO 70 
TD = FLOAT(M+1-MTCON)*DTS 
GO TO 70 


SHARP EDGE GUST (UGUST AND/OR VGUST) NOT PROVEN DUE TO 
INCONSISTENT ASSUMPTIONS REQUIRING ROTATIONAL FLOW 


XGF = T 

DO 113 L = 1,NAIRFO 
LIG = (L-1)*NODTOT 
KIG = (iil) Nea 


DO 110 IG = 1,NODTOT 

UGC IG+LIG) = 0.0 

VGC IG+LIG) = 0.0 

XG = XC IG+KIG) 

XGP1 = X( IG+KIG+1) 

IF (IG {LES NLOWERS De COmpe .120 
IF ((XGE LEP AxG CO oer 

IF CAGE ~CGE.XGr1)) COnmea i 


FAC = (AGE =" XG) / cei =) xe) 
UGC IG+LIG) = UGUST*FAC 

VGC IG+LIG) = VGUST*FAC 

GO TO 110 

UGC IG+LIG) = UGUST 

VGC IG+LIG) = VGUST 

GOFTO 110 


1P CAGE LE. Ge) Coe ioma Ye 
TP PCAGE CE 4G) COsrO 12H 


FAC = (XGF - XGP1)/(@GEeece”™ 
UG(IG+LIG) = USUST*FAC 

VG(IG+LIG) = VGUST**FAC 

GO! TO iio 

UG(IG+LIG) = UGUST 

VG(IG+LIG) = VGUST 

CONTINUE 


96 


gp (ay a a 


C762 ©) 


aan 


his 


11 


131 


lee 


141 


142 


os 


70 


CONTINUE 

IF (XGF .LE. COSALF(L)) MGUST = M 

Piers vs SNE. 0.0) GO TO 70 

IF (XGF .GT. COSALF(L)) TD = FLOAT(M+1-MGUST)*DTS 


COe1O7'C 
TRANSLATION HARMONIC OSCILLATION 


PEaeeALP . NE.w0e0)e.GO. T0012 
DO 131 I = 1,NAIRFO 


HX( 1) = DELHX(1) * SINC FREQ*T + PHA(I)) 
HYG) = DEEHY CP), s IN( FREG@4m) 

DHX(I1) = HX(I) - HXO(I) 

DHY(I) = HY(I) - HYOCI) 

UX(T) = DELHX(1)*FREQ*COS( FREQ*T+PHA( I) ) 
UY( I) = DELHY( 1 )*FREQ*COS( FREQ*T) 
CONTINUE 

XPRM = DHX(2) - DHX(1) 

PRY = DHY(2) - DHY(1) 

CALL NEWPOS(0) 

GO TO 7/0 


ROTATIONAL HARMONIC OSCILLATION 


DAL = DALP*SIN(FREQ*T) 

OMEGA( 1) = - (DALP*PI/180.) * FREQ * COS(FREQ*T) 
PESCIPHASE . EQ. 0) GO TO 14] 

ALPHA( 1)= ALPI(1) + DAL 

ALPHA(2)= ALPI(2) + DAL 

OMEGA( 2) = OMEGA(1) 

GO TO 142 

ALPHA(1)= ALPI(1) + DAL 

ALPHA(2)= ALPI(2) - DAL 


OMEGA(2) = - OMEGA(1) 

BO 143 I = 1,2 

COSALF(I) = COS(ALPHA(1)*PI/180. ) 
SINALF(I) = SIN(ALPHA(I)*PI/180. ) 
DA = ALPHA(I) - ALP(I) 
COSDA(I) = COS(DA*PI/180. ) 
SINDA(I) = SIN(DA*PI/180. ) 

DHX(1) = PIVOT * (1. -COSDA(I)) 
DHY(1I) = - PIVOT * SINDA(I) 
UY(I) = PIVOT * OMEGA(I) 
CONTINUE 


CALL NEWPOS(0) 
TRANSFORM CORE VORTEX COORDINATES W. R. T. NEW AIRFOIL POSITION 


Gero. i) GOTO 80 
DO 85 L = 1,NAIRFO 


DO 90 I=1,M-1 

XCI(L,1) = XXC(L,I) + CVVX(L,I1) * DT 

YCI(L,I1) = YYC(L,I) + CVVY(L,I) * DT 

XCO = XCI(L,1) 

YCO = YCI(L,I) 

XCI(L,I) = XCO*COSDA(L) - YCO*SINDA(L) + DHX(L) 
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Sp 1 Ie 


90 YCI(L,I) = XCO*SINDA(L) + YCO*COSDA(L) + DHY(L) 

85 CONTINUE 
CALL NEWPOS(3) 

80 CONTINUE 
WRITE (6,1001) T,DT 

1001 FORMAT (/////,' TIME STEP TK = ',F10.6,10X,'TK - TKM1 = ',F10.6,/) 
WRITE (6,1004) ALPHA( 1) ,ALPHA( 2) ,OMEGA( 1) ,OMEGA( 2) ,UX(1),UX(2), 


+ CNG) Ue2) 
1004 FORMAT (/,' ALPHA(1) = ',F10.6,5X,' ALPHA(2) = ',F10.6,/, 
+' OMEGA(1) = ',F10.6,5X, OMEGA(2) = ',F10.6, 
wena UCL) = ',F10.6,5X,'  U(2) =, 
+F10.6,/,' V(1) = " ,FiOv6 , ope ViC2) =— . Filgmee///- 
+ 1X,' NITR  VXW(1) VYW(1)  WAKE(1) THETA(1) GAMK( 1) 
+ VXW( 2) VYW(2) WAKE(2)  THETA(2) GAMK(2) ',/) 
CALCULATE THE TRAILING EDGE WAKE ELEMENT 
NUM =a) 
NITR = 0 
10° SesiseL = 1,NAIRFO 
KI = (L-1)*(NODTOT+1) 


WAKE(L) = SQRT(VYW( L)**VYW(L)+VXW(L)*VXW(L) )*DT 
THENP1(L) = ATAN2( VYW(L) ,VXW(L)) 
COSTHL(NP1+KI) = COS(THENP1(L) ) 
15 SINTHL(NP1+KI) = SIN(THENP1(L)) 
WRITE (6,1002) NITR,VXW(1),VYW(1),WAKE(1),THENP1(1),GAMK(1), 
+VXW( 2), VYW( 2) ,WAKE( 2) , THENP1( 2) , GAMK( 2) 
1002 FORMAT (15,4F10.6,E14. 6,4F10. 6,E14. 6) 
XI(NP2) = XI(NP1) + WAKE(1)*COSTHL(NP1) 
YI(NP2) = YI(NP1) + WAKE(1)*SINTHL(NP1) 
XI(NP2+1) = XI(NP1) + WAKE(2)**COSTHL( 2*NP1) 
YI(NP2+1) = YI(NP1) + WAKE(2)*SINTHL(2*NP1) 
CALL NEWPOS(1) 
CALL INFL (NITR) 
CALL COEF (NITR) 
CALL GAUSS(3,M,NITR) 


WRITE (6,*) 'ACI,NP)',(A(I,101),1=1,100) 
WRITE (6,%) 'ACI,NTOT)',(ACI,102),I=1,100) 
WRITE (6,*)'BEFORE KUTTA' 
CALL KUTTA(NITR, PVTAG) 
IF (PVTAG .LT. 0.01) GOTO 13 
NUM = NUM + 1 
WRITE (6,*)'NUM =',NUM 
IF (NUM .GT. 1) GOTO 13 
toe 1.0 
VXW(I) = 1 
7 VYW(I) =0 

connie 

13. CALL TEWAK 
DO 24 L = 1,NAIRFO 
TOL1(L) = ABS(VYW(L) - VYWK(L))/VYWK(L) 
TOL2(L) = ABS(VXW(L) - VXWK(L))/VXWK(L) 
VYW(L) = VYWK(L) 

24 VXW(L) = VXWK(L) 


IF ((TOL1(1) - Ll TOL) AND, (lOL2@i) pen om 
+. AND. (TOL1(2) . LESMGOR) AND. (ROn2C2) 0 LT, TOR) cose s20 


98 


C) Chae 
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Peeeitik .GT.. 25) THEN 
TOL = TOL*10 
MITE (6,1023) TOL 
ENDIF 
TeeenwiTR ..GT. 50) STOP 
NITR = NER e+ 1 
GO TO 10 
200) @IE.(6,1011) NITR 
1011 FORMAT (//,' CONVERGED SOLUTION OBTAINED AFTER NITR 
1023 FORMAT (//,' 7%’ TOLERANCE CRITERIA CHANGED : TOL 
CALL PRESS 
IF ((UGUST .EQ. 0.0) . AND. (VGUST .EQ. 0.0)) GO TO 300 
SIN] = SINANG(1) 


a) 
' |F10. 6) 


SIN2 = SINANG( 2) 
COS1 = COSANG( 1) 
COS2 = COSANG( 2) 
CALL FANDM(CSIN1,SIN2,COS1,COS2) 
GO TO 400 
300 SIN1 = SINALF(1) 
SINZ = SINALF( 2) 
COS1 = COSALF(1) 
COS2 = COSALF( 2) 


CALL FANDM(SIN1,SIN2,C0S1,C0S2) 
400 CONTINUE 


PEgeel TIME STEP (TADJ .NE. 0.0) IF NECESSARY 


me CADIS EQ: 020)°CO TO 95 
Were (5,2001) 
2001 FORMAT (//,' DO YOU WANT TO ADJUST TIME STEP ? 0 - NO, 1 - YES') 
PEA. (5,7) IDT 
Peerol B®, 0) GO TO 95 
DT TAdI es DT 
T TOED + DT 
WRITE (6,1006) 
1006 FORMAT (//,' BACK-TRACK COMPUTATION AND ADJUST TIME-STEP’ ,//) 


WAKE ELEMENT LEAVES TRAILING EDGE AS A CORE-VORTEX 

o5 Dae96°0L = 1;,NAaITRKFO 

CVCL,M) = SS(L)*(GAMMA(L) -GAMK(L) ) 

CVVX(L,M) = VXWC(L) 
96 CVVY(L,M) = VYW(L) 

ACIC1,M) = XICNP1) + 0. S*WAKE( 1)*COSTHL(NP1) 

TGC Use) = YIC(NP1) + 0. S*WAKE( 1)*SINTHL(NP1) 

ACTC2Z3M) = XI(NP1) + 0. S*WAKE(2)*COSTHL(NP1*2) 

TELC2..) = YICNP1) + 0. S*WAKE(2)*SINTHL(NP1*2) 


CALL NEWPOS(2) 
WRITE (6,1052) 
1052 FORMAT (//,' TRAILING VORTICES DATA',//, 
eee oe X1T CMY SX, YIIC(M)' ,6X, X1CM)' ,6X, 'Y1CM)', 
Teo CIRG! 460%, X%21(M) ,5X, YZ1€M)°, 
= 5 X2@n) sore x20M) ,e6X, CIRC2 ,/) 
DO-900 T = 1M 
Geo) Mene (6.1050) 1,XC1( 1,1), YCI(1,1),XC(1,1),Y6@1,1),CVC(1,1), 
ee? ln) YC 10 2,1), 4C(2,1),YC(2,1),6V(2,1) 


99 


1050 FORMARC1>, 1OP Pie 
CALL CORVOR 


C 
C RE-INITIALISE PARAMETERS FOR NEXT TIME STEP CALCULATION 
C 


DO 30 L = 1,NAIRFO 
Ll ( L-1)*NODTOT 


HXO(L) = HX(L) 
HYO(L) = HY(CL) 
GAMMA(L)= GAMK(L) 
ALP(L) = ALPHA(L) 
DO 30 YT = 1,NODTOT 


PHICI+L1I) = PHIKCI+L1) 
30 CONTINUE 


TOLD = T 

DT = TD 

1) = T + TD 
GO TO 40 


200 WRITE (6,1024) TOL 
1024 FORMAT (1X, /‘'%**%"* TOLERANCE CRITERION USED: TOL ' ,F10. 6) 
STOP 
END 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


C C 
C SUBROUTINE INDATA C 
C C 
C SET PARAMETERS OF BODY SHAPE C 
C FLOW SITUATION, AND NODE DISTRIBUTION C 
C C 
C USER MUST INPUT C 
C NLOWER = NUMBER OF NODES ON LOWER SURFACE C 
C NUPPER = NUMBER OF NODES ON UPPER SURFACE C 
C PLUS DATA ON BODY AND SUBROUTINE BODY C 
C C 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C 

C NACAD(I) .... NACA NUMBERS FOR THE TWO AIRFOILS 

C TAUD(I) .... MAX THICKNESS FOR THE TWO AIRFOILS 

C EPSMAD(I) ... MAX CAMBER FOR THE TWO AIRFOILS 

C PTMAXD(I) ... CHORDWISE POSN FOR MAX CAMBER FOR THE TWO AIRFOILS 

C 


SUBROUTINE INDATA 

DIMENSION TITLE( 20) 

COMMON /BOD/ IFLAG,NLOWER,NUPPER ,NODTOT, X( 202) ,Y(202), 
COSTHE( 201) ,SINTHE(201),SS(2),NP1,NP2,NP3,NP4, 
NP5 ,XSHIFT, YSHIFT,NAIRFO, X1I( 202) ,Y1(202), 
COSTHL( 201) , SINTHL( 201) 

COMMON /PARD/ NACAD(2),TAUD(2) ,EPSMAD( 2) ,PTMAXD( 2) 

READ (1,501) ITITLE 

WRITE (6,501) ITITLE 

DO 10 91 = 1:ITITLE 

READ (1,502) TITLE 

10 WRITE (6,503) TITLE 

501 FORMAT(315) 

502 FORMAT(20A4) 

503 FORMAT(1X,20A4) 


+++ 


100 


5904 FORMAT(15,2F10. 6) 
READ (1,501) IFLAG,NLOWER,NUPPER 
WRITE (6,501) IFLAG,NLOWER ,NUPPER 
READ (1,504) NAIRFO,XSHIFT,YSHIFT 
WRITE (6,504) NAIRFO,XHIFT,YSHIFT 
Bee PLAG “NE. 0) RETURN 


COMPUTATION FOR NACA 4-DIGITS SERIES 


moo 


READ (1,501) NACAD(1) 

WRITE (6,501) NACAD(1) 

ReaD (1,501) NACAD(2) 

WRITE (6,501) NACAD(2) 

DO 100 I=1,NAIRFO 

Tero NACAD(1)/1000 

IPTMAX NACAD(I)/100 - 10*IEPS 
ITAU NACAD(I) - 1000*IEPS - 100*IPTMAX 
EPSMAD( 1) Pees 0201 

PTMAXD(I) = IPTMAX*0. 1 

TAUD( I) = ITAU*0. 01 

IF (IEPS .LT. 10) GOTO 100 


COMPUTATION FOR NACA 5-DIGITS SERIES NOTING THAT TAUDC(I) IS 
COME VED AS PER 4-DIGITS SERIES. 


ep lil Dp Ls ee i ap 


PIMPADC I) = 0.2025 
Beeeao(l) = 2.6595*PIMAXD(1)**3 
100 CONTINUE 
RETURN 
END 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


BODY AND NACA 


C C 
C SUEBREUITRE SETUP C 
C C 
C SETUP COORDINATES OF PANEL NODES AND SLOPES OF PANELS C 
C COORDINATES ARE READ FROM INPUT DATA FILE UNLESS C 
C THE AIRFOIL IS OF NACA XXXX OR NACA 230XX TYPE C 
C C 
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
C 

C NACACT) .... DUMMY NACA NUMBER FOR TRANSFER TO SUBROUTINES 

C BODY AND NACA 

C TAUC I) .... DUMMY MAX THICKNESS FOR TRANSFER TO SUBROUTINES 

C BODY AND NACA 

C EPSMAX(I) .... DUMMY MAX CAMBER FOR TRANSFER TO SUBROUTINES 

C BODY AND NACA 

C PTMAX(I) .... DUMMY CHORDWISE POSN FOR TRANSFER TO SUBROUTINES 

C 

C 


SUBROUTINE SETUP 

COMMON /BOD/ IFLAG,NLOWER ,NUPPER ,NODTOT,X(202) ,Y¥(202), 
COSTHE( 201) ,SINTHE(201) ,$S(2) ,NP1,NP2,NP3,NP4, 
NP5 , XSHIFT, YSHIFT ,NAIRFO,XI(202) ,YI(202), 
COSTHL( 201) ,SINTHL( 201) 

COMMON /PARD/ NACAD(2),TAUD(2) ,EPSMAD(2) , PTMAXD( 2) 

COMMON /PAR/ NACA,TAU,EPSMAX, PTMAX 


+++ 


101 


C2. GaGa 


Le 1 Dal Sp 


Cat) 


i C2) <2 


100 


110 


120 


10 


pol 
210 


COMMON /NUM/ PI,PI2INV 
Ea = 3.1415926585 
PIZINV = 75/7 FP. 


SET COORDINATES OF NODES ON BODY SURFACE 


DO 210 I=1,NAIRFO 
NODTOT = NLOWER + NUPPER 
K=(I-1)*(NODTOT+1) 

IF (IFLAG .NE. 0) GO TO 10 


NACA SERIES AIRFOIL CALCULATIONS 


NACA=NACAD(T) 
TAU=TAUD(T) 
EPSMAX=EPSMAD( TI) 
PTMAX=PTMAXD(T) 


NPOINT = NLOWER 

SIGN = -1.0 

NSTART = 0 

BOeLTIO NSURE = si2 

BOri00 N= 1,NPOIND 

FRACT = FLOAT(N-1)/FLOATCNPOINT) 
Z =“.o°Cl. - COsCPItrRaGm,) 
J = NSTART + N 


CALL BODY(Z,SIGN,XD,YD) 
XI(J+K) = XD 

YICJ+K) = YD 

CONTINUE 
NPOINT 
SIGN 
NSTART 
CONTINUE 
XIC I*(NODTOT+1) ) 
YiC1*(NOLTOT AIS 
GOv10 210 


NUPPER 
120 
NLOWER 


Crh) 
YIGisk) 


AIRFOIL DATA INPUT THAT ARE NOT NACA SERIES AIRFOIL. 


READ (1,501) (XI(J+K) ,J=1,NODTOT+1) 
READ (1,501) (YI(J+K) ,J=1,NODTOTH1) 
WRITE (6,501) (XI(J+K) , J=1,NODTOT+1) 
WRITE (6,501) (YIC(J+K) , J=1,NODTOT+1) 
FORMAT (6F10. 6) 


CONTINUE 
NPT = NODTOT + 1 
NP2 = NAIRFO*NP1+1 


SET SLOPES OF PANELS AND CALCULATE AIRFOIL PERIMETER 


DO 220 I = 1,NAIRFO 

K = (I-1)*(NODTOT+1) 
SS(I) = 0.0 

DO 200 J = 1,NODTOT 

DX = XI(J+1+K) - XI(J+K) 
DY = YI(J+1+K) - YI(J+K) 


200 
220 


DIST SQRT(DX*DX +DY*DY) 
Sou L) SO ee ILL Se 
SINTHL( J+K) DY/DIST 
COSTHL( I+K ) DX/ DIST 
CONTINUE 

CONTINUE 

Fe tURN 

END 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCcccccccccccccccccccceccccccccccccce 


C 
C 
C 
C 
C 
C 
C 
C 
C 
C 


SUBROUTINE BODY(Z,SIGN,X,Y) 


RETURN COORDINATES OF POINT ON THE BODY SURFACE 
Z= 
X,Y = CARTESIAN COORDINATES 
SIG 


NW = +1. FOR UPPER SURFACE 
-1. FOR LOWER SURFACE 


C 

C 

C 

C 

C 

NODE-SPACING PARAMETER C 

; C 
I C 
C 

C 

C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


SUBROUTINE BODY(Z,SIGN,XD, YD) 
COMMON /PAR/ NACA, TAU, EPSMAX, PTMAX 
PP (SIGMes LT. 0.0) Lo Wee 
CALL NACA45(Z,THICK, CAMBER, BETA) 


XD = Z - SIGN*THICK*SIN( BETA) 
YD = CAMBER + SIGN*THICK*COS( BETA) 
RETURN 
END 
GemeeeecoeccCCCCCCCcceeeeccCCCECECCCCCCCCCCECCeCcececececcccccecccccccccccc 
c 6 
@ SUBROUTINE NACA4S(Z, THICK, CAMBER, BETA) C 
E e 
@ EVALUATE THICKNESS AND CAMBER € 
C FOR NACA 4- OR 5-DIGIT AIRFOIL C 
c E 
pao eee esos eee CCCCECCCCCCECCCCCCCCCECCCCCCCEeccEecgeececcccceeccceccccccccc 
° 
@ Te ABSCISSA OF CAMBERLINE POINT 
¢ ae LOCAL CARTESIAN COORDINATES OF NACA-AIRFOIL 
6 SIGW ...... SURFACE INDICATOR: SIGN = 1, UPPER SURFACE 
& =-1, LOWER SURFACE 
E THC .. 2... THICKNESS AT THE POSITION Z 
c CAMBER ...... CAMBER AT THE POSITION Z 
C Pom ...... ANGLE BETWEEN CAMBER LINE AND X-AXIS AT POSITION Z 
C fy MAXIMUM THICKNESS (INPUT) 
C EPSMAX ...... MAXIMUM CAMBER (INPUT) 
C PTMAX ...... COORDINATE POSITION OF MAXIMUM CAMBER ( INPUT) 
E 
SUBROUTINE NACA4S5(Z, THICK, CAMBER, BETA) 
COMMON /PAR/ NACA,TAU,EPSMAX, PTMAX 
THICK = 0.0 
TZ TMI, Do P= 110) lomo aor, 
THICK = 5.*TAU*(.2969*SQRT(Z) - Z*(.126 + Z*(. 3537 
+ Ae een jt) 
100 IF (EPSMAX .EQ. 0.0) GO TO 130 


103 


qa CC) ©) 


lao 


110 
120 
130 


Oaqn 


140 


qaqa 


qan 


Beye) 


IF (NACA .GI.> 9993)3G0 TO 140 
CAMBERLINE OF NACA 4-DIGIT SERIES 
IF (2. Glas Pierce tCerre 
FORWARD PART OF CAMBER LINE 


CAMBER = EPSMAX/PTMAX/PTMAX*(2. *PTMAX - Z)*Z 
DCAMDX = 2. *EPSMAX/PTMAX/PTMAX*(PTMAX - Z) 
GOs LO. 120 


AFT PART OF CAMBERLINE 


CAMBER = EPSMAX/(1. -PTMAX)**2*(1. + Z - 2. *PTMAX)*(1. - Z) 
DCAMDX = 2. *EPSMAX/( 1. -PTMAX)**2*(PTMAX - Z) 

BETA = ATAN(DCAMDX) 

RETURN 

CAMBER = 0.0 

BETA = 0.0 

RETURN 


CAMBERLINE OF NACA 5-DIGIT SERIES 
IE (2 .Gi2Pnwis® GO TO 150 
FORWARD PART OF CAMBER LINE 


W = Z/PTMAX 

CAMBER = EPSMAX*W*((W - 3. )*W + 3. - PTMAX) 
DCAMDX = EPSMAXK*3.*W(1. - W)/PTMAX 

GOTO 7120 


AFT PART OF CAMBERLINE 


CAMBER = EPSMAX*(1. - Z) 
DCAMDX = - EPSMAX 

GO TO 120 

END 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 


SUBROUTINE GAUSS(CNRHS ,M,NITR) 


SOLUTION OF LINEAR ALGEBRAIC SYSTEM BY 
GAUSS ELIMINATION WITHOUT PARTIAL PIVOTING 


°A = COEFFICIENT MATRIX 
NEQNS = NUMBER OF EQUATIONS 
NRHS = NUMBER OF RIGHT HAND SIDES 


RIGHT-HAND SIDES AND SOLUTIONS STORED IN 
COLUMNS NEQNS+1 THRU NEQNS+NRHS OF °A 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


SUBROUTINE GAUSS(NRHS ,M,NITR) 
COMMON /BOD/ IFLAG,NLOWER ,NUPPER,NODTOT, X( 202) ,Y(202), 


104 


C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 


COSTHE(201) ,SINTHE( 201) ,SS(2) ,NP1,NP2,NP3,NP4, 
NES oP ey OHIPT | NAIRFOsX1G202),Y10202).. 
COSTHL( 201) ,SINTHL( 201) 

COMMON /COF/ AC201,211),KEQNS 


+++ 


° WRITE (9,900) M,NITR 
©9708) FORMAT (1X, °M ='/13,' NITR =' ,13/) 
C Wel te"C9,910)(1,AC1,NP3) ,ACI1,NP4) ,AC1,NP5) , 1=1, 2*NODTOT) 


OO FORMA (1X ,15,3E14. 6) 
IF (M.EQ.0) KEQNS = NODTOT 
NEQNS KEQNS * NAIRFO 
NP NEQNS + 1 
NTOT NEQNS + NRHS 
re CN pol. 0), COMER GO 


ei tt ott 


GAUSS REDUCTION 


aang 


BO150» I 
IM = 


ELIMINATE (I-1)TH UNKNOWN FROM 
ITH THRU (NEQNS)TH EQUATIONS 


C27 Cr OQ 


wo 150° J = 1,NEQNS 
R AC J,1IM)/ACIM, IM) 
DO 150 =" aNor 

8 ACI,K) bGiGh) ~ RMS CIM;K) 
Sa TO 170 


Wa tl 


GAUSSIAN ELIMINATION ON ONLY THE RIGHT-HAND-SIDES 


Oanqn 


160 DO 180 

eM 

DO 180 

R 

DO 180 
730 GacJ,K) 
170 CONTINUE 


2 ,NEQNS 

- 1 

I ,NEQNS 
A(J,IM)/ACIM, 1M) 

= NP NTOT 

A(J,K) - R*ACIM,K) 


unit GQ ite 
We tl 


BACK SUBSTITUTION 


Oqan 


DO 220 K = NP,NTOT 
A(NEQNS ,K) = A(NEQNS,K)/A(NEQNS , NEQNS ) 
DO 210 L = 2,NEQNS 
I NEQNS + 1 - L 
IP ieee 
DO 200 J = IP,NEQNS 
200 A(I,K) ACI,K) - ACI,J)*A(J,K) 
210 GT K) A(1I,K)/A(I,1) 
220 CONTINUE 
RETURN 
END 
EiGeG eee serceGeGEeGGeCceEscGCCeCCCCCCCCECCEGCCCCGGEGEGECGCCCCCCEEECCCCCCCC 
C é 
SUBROUTINE VELDIS(SINALF ,COSALF) 


Hueah tt 


C C 
C C 
C COMPUTE STEADY FLOW PRESSURE DISTRIBUTION AND VELOCITY C 


105 


C POTENTIAL AT MID-POINTS OF PANELS FOR THE STEADY FLOW CASE C 
C C 
cccecccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
C 

C VIANG. 222228 TANGENTIAL VELO COMP OF A FIXED POINT OF THE MOVING 

C LOCAL FRAME OF REFERENCE. 

C PHL lpia oe DISTURBANCE VELO POTENTIAL AT THE MID POINT OF THE 

C I-TH PANEL 

C PHILE(L).... DIFFERENCE OF THE POTENTIALS OF THE LEADING EDGE TO 

C THE LOWER TRAILING EDGE FOR THE RESPECTIVE AIRFOIL 

C PIN(L) .... DIFFERENCE OF THE POTENTIALS AT A POINT 1000 CHORD 

C LENGTH UPSTREAM OF THE LE FOR THE RESPECTIVE AIRFOIL 

C SUMC(L) .... GAMMA ASSOCIATED WITH THE INTEGRATION OF THE DISTURB- 
C ANCE VELOCITY AROUND THE WHOLE AIRFOIL 

C 


SUBROUTINE VELDIS 
COMMON /BOD/ IFLAG,NLOWER,NUPPER ,NODTOT, X(202),Y(202), 
COSTHE( 201) ,SINTHE( 201) ,SS(2),NP1,NP2,NP3,NP4, 
NP5 ,XSHIFT, YSHIFT ,NAIRFO, XI(202) ,YI(202), 
COSTHL( 201) , SINTHL( 201) 
COMMON /COF/ A(201,211),KEQNS 
COMMON /CPD/ CP(200),SCL,T,SCM,SGAM 
COMMON /NUM/ PI,PI2INV 
COMMON /SING/ Q(200) ,GAMMA(2),QK( 200) ,GAMK(2) 
COMMON /POT/ PHI( 200) ,PHIK( 200) 
COMMON /GUST/ UG(200) , VG( 200) ,XGF ,UGUST, VGUST 
COMMON /EXTV/ UE(200),VN(200) 
COMMON /GEOM/ SINALF(2),COSALF( 2) ,OMEGA(2),UX(2),UY(2),PIVOT, 
A XPRM , YPRM 
DIMENSION CKUTTA( 2), PHITEL(2) , PHITEU(2) 
7770 DIMENSION CONTR2(2) 
6666 REAL *8 PIN(2),VELX 
DIMENSION WGHT(5),PLOC(5),SUMC(2), 
+AANP1(50,50,5),AANP2(50,50,5),BBNP1(50,50,5),BBNP2(50,50,5), 
+COSTHP( 102,6),SINTHP( 102,6) , AANP4( 2) , PHIL( 102) ,UGU( 100,6), 
+VGU( 100 ,6) , PHILE(2) 


+++ 


C 
C LOCATION AND WEIGHTING VALVES FOR THE GAUSSIAN QUADRATURE USED 
C TO INTEGRATE FOR THE VELO POTENTIAL AROUND THE AIRFOIL 
C 

DATA WGHT/. 11846344,.23931434,. 28444444, 

+ » 23931434,.11846344/ 
Pei PLOC/. 04691008,. 23076535,. 50000000, 
; /69234606,, Jaa0Uoog, 

C 
C FIND VT AND CP AT MID-POINT OF I-TH PANEL 
C 

DO 140 L= 1,NAIRFO 

KI = (L-1)*NP1 

LI = (L-1)*NODTOT 

SUMC(L) = 0.0 

DO 130 I = 1,NODTOT 

AMID = .5*(X(I+KI) +X(I+KI+1)) 

YMID = ,5*(YCI+KI]) +YC4Ki Ge 

DX = XC IFKIF1) XG 

DY = Y(I+KI+)l)-YCake 


106 


CC? Cl Ci C1 Gre @) 


OrIOaAIANN 


100 


120 
P50 


10 
140 


mast SQRT( DX*DX+DY*DY ) 


VIANG = COSALF(L)*COSTHL(I+KI) + SINALF(L)*SINTHL(I+KI) 
VNORM = SINALFC(L)*COSTHL(CI+KI) - COSALF(L)*SINTHL(I+KI) 
WIFREE = VTANG 
VACT = VTANG 


ADD CONTRIBUTION OF J-TH PANEL 


1. CONTRIBUTION OF J-TH PANEL TO THE VELO COMP OF THE MIDPT 
OF TH I-TH PANEL 


2. CONTRIBUTION TO THE VELO POTENTIAL. THIS IS DONE BY 
INTEGRATING OVER SMALLER PANELS OF THE AIRFOIL. 


BO 155 F= 1,5 

DX = PLOC(K)*( XC I+KI4+1)-X(CI+KI)) 
DY = PLOC(K)*( YC I+KI+1)-YC1I+KI)) 
XMID = X(I+KI) + DX 

YMID = Y(I+KI) + DY 

VDUM = 0.0 

DO 150 LM = 1,NAIRFO 

KJ = (LM-1)*NP1 

LJ = (LN-1)*NODTOT 

DO 120 J = 1,NODTOT 

FLOG = 0.0 

FTAN = PI 


tEeCorRS .EQ. IKI) GO TO 100 


Dad = hi1D - AC JtKJ) 

DXJP =p e- A J2RI+1) 

pe J ="SHID - X(JtKI) 

DYruP = ERID - 3 J-3J41) 

FLOG = . MAEOG( (DAIP“DEIPTDYIPADY IP ) /( DXJ*DXJ+DYJ*DYJ) ) 

FTAN = 2 eee ( Be JP Dx J -DxgP4Dr J , DAJP*DxXJ+DY JP*DYJ) 

GHP i = COST -( IF] )-COSTHE( J+KJ) + SINTHECI+K1)*SINTHE( J+KJ) 
Sete «= SHNWRHE( I+KE )*COSTHE( J+KJ) - COSTHE(C I+KI)*SINTHE( J+KJ) 
AA = PR2ZING+(FPAN'CTINTS + FLOG*STIMTJ) 

B = PIZINV*(FLOG*CTINTJ - FTAN*STIMTJ) 

VDUM = VoUM - B*OQ(J+LJ) +GAMMA(LM)*<AA 


TR(K .EQ. 3) VACT = VACT = B*Q(J+LJ) + GAMMA(LM)*AA 
IM K .EQ. 3) VNORM = WRORM + AA*Q(J+LJ) + GAMPIACLH)=8B 
CONTINUE 


CONTINUE 

VIANG = VIANG +VDUM *WGHT(K) 
SUMC(L) = SUMC(CL) +VDUM*DIST*WGHT(K) 
CONTINUE 


PHICI+L1) = (VIANG-VIFREE)*DIST 
GP(T+b1)= 1.0 - MACT“VACT 

UEC I+LI)= VACT 

VNC I+L1I)= VNORM 

CONTINUE 

CONTINUE 


COMPUTE DISTURBANCE POTENTIAL AT THE LEADING EDGE BY LINE 
INPEGRAL @F THE VELOCITY #1ELD 
FROM UPSTREAM (AT INFINITY) TO THE LEADING EDGE 


po 55 L = 1,NAIRFO 
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202 © 


qa 


i 1 Ge 


20 
40 


30 
55) 


230 


YMID 
XLE 
XL 

XL 
NPHI 
PIN(L) 
DO 30 
FRACT 
XLP 

LP Gie.E 
XLP 
DELX 
XMID 
XMID 
Paulie, 
XL 
VELX 


nnn ntthonrtiei t tien 


ADD CONTRI 


DO 40 ~LN 
KJ 

a 

DO 20 
DXJ 

DXJP 

DYJ 

DYJP 
FLOG 
FTAN 
CALMTJ 
SALMTJ 
APY 

BPY 

VELX 
CONTINUE 
CONTINUE 
PIN(L) = 
CONTINUE 
CONTINUE 


munud nnn nnd nah ll 


COMPUTATIO 


DO 240 L 
LI 
Por 


BEGIN WITH 


DOeZ30) 1 
PHC 
Farell) 
PHP 
PHITEL(L) 


PIVOT * SINALFCL) + (L=1)4¥shiee 
-PIVOT*COSALF(L)+( L-1)*XSHIFT 
XLE 
0.0 
10 * NLOWER 
070 
= 1,NPHI 
FLOAT(1I)/FLOAT(NPHI ) 
-100.0 * (1.0 - COS(0. 5*PI*FRACIRY -xie 
1) XLP = -0. 000197+XLE 
-10.0.*% (1 Os="C0S (G7 5*PI*ERACT)D 
XL =3xEP 
O. 5°*(XL+XLP) 
O. 5**(XL+XLP)*COSALF(L) 
O. 5*(XL+XLP)*SINALF(L) 
XLP 
UGUST 


BUTION OF J-TH PANEL 


= 1,NAIRFO 

(LN-1)*NP1 

( LN-1)*NODTOT 

= 1,NODTOT 

XMID - X(J+KJ) 

XMID - X(J+KJ+1) 

YMID - Y(J+KJ) 

YMID - Y(J+KJ+1) 

. 5*ALOG( (DXJP*DXJP+DYJP*DY JP) / (DXJ*DXJ+DYJI*DYJ) ) 
ATAN2( DYJP*DXJ-DXJP*DYJ , DXJP*DXJ+DYJP*DYJ) 
-COSALF(L)*COSTHE( J+KJ) - SINALF(L)*SINTHE( J+KJ) 
-SINALF(L)**COSTHE( J+KJ) + COSALF(L)*SINTHE(J+KJ) 
PI2INV*(FTAN*CALMTJ + FLOG*SALMTJ) 
PI2INV**(FLOG*CALMTJ - FTAN*SALMTJ) 

VELX - DPROD(BPY,Q(J+LJ)) +DPROD(GAMMAC( LN) , APY) 


PIN(CL) + VELX * DBLE(DELX) 


N OF THE VELOCITY POTENTIAL FOR MIDPOINT OF EACH PANEL 


1,NAIRFO 
(L-1)*NODTOT 
-PIN(L) 


Hou tl 


LOWER SURFACE 


= NLOWER,1,-1 
PHP-PHI(I+LI) 
0. 5*(PHP+PHC) 
PHC 
PHC 


Hou uu 
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C RESET FOR UPPER SURFACE 


C 
PHP = -PIN(L) 
Peecoo «| = NLOWER+1 ,NODTOT 
PHC = PHP+PHI(I+LI) 
PHICI+LI) = 0. 5*( PHP+PHC) 
z50 PHP = PHC 
Peirevocn) = PHC 


240 CONTINUE 
DO 334 L =1 , NAIRFO 
WRITE (6,1010) L 
WRITE (6,1000) 


KI = (b-1) 7 "Net 

LI = (L-1) * NODTOT 

SUMC(L) = SUMC(L)/SS(L) 

DO 333 I = 1 , NODTOT 

XMID = 2 5 * CXC a Gti +1) 
ELD Sono YC) te eR) ) 


orn oO, 1050) ItL1SXMID, YMID, Q(1+L1) ,GAMMA(L) ,€P(1+L1) ,UECI+L1), 
+ VNC I+L1),PHIC I+L1) , SUMC(L) 

334 WRITE (6,235) PIN(L) 

235 FORMAT (1X,'PHI AT LEADING EDGE = ',F10. € 6,/) 

1O@0 FORMAT(/,4X,'J' ,4X,'XCJ)" 5x ‘YCJ) lone 'O(J)' 6X, 'GAMMA' ,5X, 
- NepCd)' ,6X,'V(J)' .6X, ¥NCJ) 56x, 'PHI' 4X." INTGAMMA' , /) 

1010 FORMAT(//' AIRFOIL NUMBER' ,14,/) 

1050 FORMAT(15,9F10. 6) 


RETURN 
END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C 

C SUBROUTINE FANDM(SINALF , COSALF ) C 

C C 

C COMPUTE AND PRINT OUT CD,CL,CM C 

C INTEGRATE PRESSURE DISTRIBUTION BY TRAPEZOIDAL RULE C 

C C 

cecceccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccce 

C 

C Cle 1) aera PRESOURE COEFFICIENT OF THE 1-TH PANEL 

C CUS aa COREEICIENT OF LIFT FOR THE L-TH AIRFOIL 

C C1 ee COEFFICIENT OF DRAG FOR THE L-TH AIRFOIL 

C GMC) ....... COEFFICIEN® OF MOMENT FOR THE L-TH AIRFOIL WITH 

C RESPECT TO THE LEADING EDGE 

C OC) eee a COEFFICEIENT OF TOTAL FORCE IN X-DIR OF GLOBAL SYS. 

C GAC) 2.0506 COEFFICEIENT OF TOTAL FORCE IN Y-DIR OF GLOBAL SYS. 

C 

C 


SUBROUTINE FANDM(SIN1,SIN2,COS1,COS2) 

COMMON /BOD/ IFLAG,NLOWER ,NUPPER ,NODTOT, X(202),Y(202), 
COSTHE( 201) ,SINTHE( 201) ,SS(2),NP1,NP2,NP3,NP4, 
NPS ,XSHIFT , YSHIFT ,NAIRFO,XI( 202) ,Y1I( 202), 
COSTHL( 201) , SINTHL( 201) 

COMMON /CORV/ CV(2,200) ,XC(2,200),YC(2, 200) ,M,TD,CVVX(2,200), 

CVVY( 2,200) ,XCI(2,200) ,YCI(2,200) 

antes /CPD/ CP(200),SCL,T,SCM, SGAM 

COMMON /SING/ Q( 200) ,GAMMA( 2) ,QK( 200) ,GAMK(2) 

COMMON /GEOM/ SINALF(2),COSALF(2),OMEGA( 2) ,UX(2) ,UY(2) , PIVOT, 


+++ 


109 


Ci 


aan 


De al a 1 ep 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


©) Qt Circa. C22 C2 C2 


XPRM, YPRM 
DIMENSION CM(2),CD(2),CL(2),CFX(2) ,CFY(2) 


INITIALISE COEFFICIENT 
DO 110 L = 1,NAIRFO 
CM(L) = 0.0 
GEXGL)  — Ome 
110 CFY(L) = 0.0 
INTEGRATE AROUND THE AIRFOIL TO GET GLOBAL X AND Y FORCES 
DO 120 L = 1,NAIRFO 
KI = (L-1)*(NODTOT+1) 
LI = (L-1)*NODTOT 
DO 100 I = 1,NODTOT 
XMID = .5%*(XI(I+KI) + XI(I+KI4#1)) 
YMID = .5*(YI(I+KI) + YI(I+KI+#1)) 
DX = XI(I+KI+1) - XICI+KI) 
DY = YI(I+KI+1) - YICI+KI) 
CFX(L) = CFX(L) + CP(I+LI)*DY 
CFY(L) = CFY(L) - CP(I+LI)*DX 
CM(L) = CM(L) + CP(I+LI)*(DX*XMID + DY*YMID) 
100 CONTINUE 
120 CONTINUE 
DECOMPOSE INTO LIFT AND DRAG COMPONENTS W.R.T. RESP LOCAL SYSTEM 
DO 130 L =1,NAIRFO 
CD(L) = CFX(L)*COSALF(L) + CFY(L)*SINALF(L) 
CL(L) = CFY(L)*COSALF(L) - CFX(L)*SINALF(L) 
130 WRITE (6,1000) b CDG eon menc.) 
IF (M .EQ. 0) RETURN 
IF (SCL NE 020) sqhnN 
WRITE (8,1100) T,CL(1)/SCL,CLC2)/ Str eM@inn SCMNG Maen SCM Chains 
+CD(2),GAMK(1)/SGAM, GAMK( 2) /SGAM 
ELSE 
WRITE €8y1100) T,CL(1),CL(2), CMG ONC ena) .cn(2) 
ENDIF 
1000 FORMAT(//,' AIRFOIL NO ',1I4,/,' CD =',F10.6, 
He Cl= SFiGs6o CM=" ,F10. 6) 
1100 FORMAT(9F10. 6) 
RETURN 
END 


SUBROUTINE INFL (NITR) 


CALCULATE INFLUENCE COEFFICIENTS 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


INFLUENCE COEFFICIENTS ON THE AIRFOIL DUE TO THE AIRFOIL: 
AAN( I,J) .... NORMAL VELO AT THE MIDPOINT OF THE I-TH PANEL DUE 


TO A SOURCE-=DIST OP UNIT SIRENGTH ON THESI-Ti PANEL 


SUMAAN(I,L) .. NORMAL VELO AT THE MIDPOINT OF THE I-TH PANEL DUE 


110 


C2). C2: G2 oC) C2 


i a SB | <p 7 > Jie ie Tk > YL > | a Ug am TA a> TR > at a a De J > Je > Yk a J > I ae Fe a J ae Tk ep Ym ge > | gp gg ie | >) 


C2 2 CaO) 


+++ 


SOURCE DIST OF UNIT STRENGTH FROM THE L AIRFOIL 
BBN( I,J) .... NORMAL VELO AT THE MIDPOINT OF THE I-TH PANEL DUE 

TO A VORTEX-DIST OF UNIT STRENGTH ON THE J-TH PANEL 
SUMAAN(I,L) .. NORMAL VELO AT THE MIDPOINT OF THE I-TH PANEL DUE 

VORTEX DIST OF UNIT STRENGTH FROM THE L AIRFOIL 


INFLUENCE COEFFICIENT ON THE WAKE ELEMENT: 

AYNP1(L,J) ... Y-VELO COMP AT THE MIDPOINT OF THE WAKE PANEL FROM 
THE L-TH AIRFOIL DUE TO A SOURCE DIST OF UNIT STRE- 
NGTH FROM THE J-TH PANEL 

AYNP1(L,NP3) . Y-VELO COMP AT THE MIDPOINT OF THE WAKE PANEL FROM 
THE L-TH PANEL DUE TO A SOURCE DIST OF UNIT STREN- 
GTH FROM THE WAKE PANEL OF THE OTHER AIRFOIL. 
(USED ONLY FOR BXNP1(L,NP3) SINCE THERE IS NO SOURC 
DIST ON THE WAKE PANEL) 

BYNP1(L,J) ... Y-VELO COMP AT THE MIDPOINT OF THE WAKE PANEL FROM 
THE L-TH AIRFOIL DUE TO A VORTEX DIST OF UNIT STRE- 
NGGUPeROM THE J-1h PANEL 

BYNP1(L,NP3) . Y-VELO COMP AT THE MIDPOINT OF THE WAKE PANEL FROM 
THE L-TH AIRFOIL DUE TO A VORTEX DIST OF UNIT STR- 
ENGTH FROM THE WAKE PANEL OF THE OTHER AIRFOIL. 

GryveiCL,N) ... Y=VELO COMP AT THE MIDPOINT OF THE WAKE PANEL FROM 
THES i - Tian POmL DUE TO THE N-TH CORE VORTEX OF 
UNIT STRENGTH 

CAePicL,N) ... X=VELO COMP AT THE MIDPOINT OF THE WAKE PANEL FROM 
THE G-TH AIRFOIL DUE TO THE N-TH CORE VORTEX OF 
UNIT STRENGTH 


INFLUENCE COEFFICIENTS ON THE AIRFOIL DUE TO THE WAKE 


SWMNCGNGI) ... NORMAL VELO AT THE MIDPOINT OF THE I-TH PANEL DUE 
TO ALL POINT VORTICES OF ACTUAL STRENGTH. 
SUMCCT(I) ... TANGENTIAL VELO AT THE MIDPOINT OF THE I-TH PANEL 


DUE TO ALL POINT VORTICES OF ACTUAL STRENGTH. 


SUBROUTINE INFL (NITR) 

COMMON /BOD/ IFLAG,NLOWER ,NUPPER ,NODTOT,X( 202) ,Y( 202), 
COSTHE(201), SINTHE(201), ssn MENEZ, NPS, NES, 
NP5 ,XSHIFT, YSHIFT ,NAIRFO, XI( 202) ,YI(202), 
COSTHL( 201) ,SINTHL( 201) 

COMMON /NUM/ PI,PI2INV 

COMMON /WAK/ VYW(2),VXW(2) ,WAKE(2) , DT 

GOMMON /CORV/ CV(2,200),XC( 2,200) ,YC(2,200) ,.M,ID,CVVX( 2,200), 


+ GVViG2 200), ACI( 2, 2000,,YCI(2, 200) 


COMMON /INFI/ AAN(201,201) ,BBN(201,201),AYNP1(2,201),BYNP1(2,201), 


+ SUMAAN( 201, 2) ,SUMBBN(201,2) 


COMMON /INF2/ SUMCCN( 201) ,SUMCCT( 201) ,CYNP1(2,200), 


+ CXNP1( 2,200) 


COMMON /GEOM/ SINALF(2), CGOSALP (2), OMEGAC2) , UX(CZ) ,UY(2), PIVOT, 


+ ABRMew PRM 


DIMENSION GAYNP1(2,201) ,GBYNP1(2,201) 
IF @ .GT. 1) GO TO 500 


INFLUENCE COEFFICIENT ON THE I-TH PANEL BY THE J-TH PANEL FROM 
THE SAME AIRFOIL. 


111 


qaan 


100 


110 
120 
200 
5119) 


DO 200 L = 1,NAIRFO 

LI = (L-1)*NODTOT 

KI = (L-1)*NP1 

DO 120 I = 1,NODTOT 

XMID =. 9° (XPCrte) + XTC 
YMID = ,5*(YI@RERT) + Yigg 


SUMAANCI+LI,L) = 0.0 
SUMBBNCI+LI,L) = 0.0 
DO 110 J = 1,NODTOT 


FLOGuN = Ome 

FTAN = PI 

IF (J+LI .EQ. I+LI) GO TO 100 

DXJ = XMID - XI(J+KI) 

DXJP = XMID - XI(J+KI+1) 

DYJ = YMID - YI(J+KI) 

DYJP = YMID - YI(J+KI+1) 

FLOG = .5*ALOG((DXJP*DXJP+DYJP*DYJP) /( DXJ*DXJ+DYJ*DYJ) ) 

FTAN = ATAN2(DYJP*DXJ-DXJP*DYJ, DXJP*DXJ+DYJP*DYJ) 

CTIMTJ = COSTHL( I+KI)*COSTHL(J+KI) + SINTHL( I+KI)*SINTHL(J+KI) 
STIMTJ = SINTHL(I+KI)*COSTHL(J+KI) - COSTHL( I+KI)*SINTHL(J+KI ) 


AANCI+LI ,J+L1) 
BBNCI+LI , J+L1) 
SUMAANCI+LI1,L) 
SUMBBNCI+LI1 ,L) 
CONTINUE 
CONTINUE 
CONTINUE 
CONTINUE 
IF (NITR {SiO COMO e 27a 


PI2ZINV*(FTAN*CTIMTJ + FLOG*STIMTJ) 
PIZINV*( FLOG*CTIMTJ - FTAN*STIMTJ) 
SUMAANCI+LI,L) + AANCI+LI ,J+LI ) 
SUMBBNC I+LI,L) + BBNCI+LI,J+LI) 


INFLUENCE COEFFICIENT ON THE I-TH PANE BY THE J-TH PANEL FROM 
THE OTHER AIRFOIL. 


DO 270 L = 1,NAIRFO 

LI = (L-1)*NODTOT 

KI = (L-1)*NP1 

DO 260 I = 1,NODTOT 

XMID = .5*(XCI+KI) + XC I+KI+1)) 
YMID = . 5 *CYCI+KI] ) +e v1 +h) 


SUMAANCI+L1I,3-L) = 0.0 
SUMBBN(I+LI,3-L) = 0.0 


DO 250 J = NODTOT+2,2*NODTOT+1 

DXJ = XMID - X(J-KI) 

DXJP = XMID - X(J-KI+41) 

DYJ = YMID - Y(J-KI) 

DYJP § = YMID - Y(J-KI+1) 

FLOG = .5*ALOG((DXJP*DXJP+DYJP*DYJP) /( DXJ*DXJ+DYJ*DYJ) ) 
FTAN = ATAN2(DYJP*DXJ-DXJP*DYJ , DXJP*DXJ+DYJP*DYJ) 

CTIMTJ = COSTHE(I+KI)*COSTHE(J-KI) + SINTHE( I+KI)*SINTHE(J-KI) 
STIMTJ = SINTHE(I+KI)*COSTHE(J-KI) - COSTHE( 1+KI)*SINTHE( J-KI) 


AANCI+LI,J-LI-1) = PI2ZINV*(FTAN*CTIMTJ + FLOG*STIMTJ) 
BENCH J - br) PI2ZINV*(FLOG*CTIMTJ - FTAN*STIMTJ) 
SUMAAN(CI+LI ,3-L) SUMAANCI+LI,3-L) + AANCI+LI,J-LI-1) 
SUMBBN(I+LI1,3-L) SUMBBN(I+L1,3-L) + BBNCI+LI,J-LI-1) 


250 CONTINUE 
260 CONTINUE 


lee 


aqany 


aqan 


gqaan 


270 


130 


300 


CONTINUE 


END OF STEADY FLOW CALAULATION FOR INFLUENCE COEFFICIENT. 


im CiPEQ.0) RETURN 
CONTINUE 


INFLUENCE COEFFICIENT ON THE WAKE ELEMENT FROM THE AIRFOIL. 


DO 130 L = 1,NAIRFO 

I = NP1*L 

XMID = .57(X(1) + X(NAIRFO*NP14L) ) 

YMID = .5°*(Y(1) + Y(NAIRFO*NP1+L)) 

DO 130 MM = 1,NAIRFO 

108) =(MM-1)*NODTOT 

ea =(MM-1)**NP1 

bowie J = 1,NODTOT 

Do = XMID - X(J+KJ) 

DXJP = XMID - XC J+KJ+1) 

DYJ = YMID - Y(J+KJ) 

DYJP = YMID - Y(J+KJ+1) 

FLOG = . 5*ALOG((DXJP*DXJP+DYJP*DYJP) /( DXJ*DXJ+DYJ*DYJ) ) 
FTAN = ATAN2(DYJP*DXJ-DXJP*DYJ , DXJP*DXJ+DYJP*DYJ) 
CTIMTJ = COSTHE(I)?*COSTHE( J+KJ) + SINTHE(1)?*SINTHE(J+KJ) 
STIMTJ = SINTHE(1)**COSTHE(J+KJ) - COSTHE(1)**SINTHE( J+KJ) 


GAYNP1(L,J+LJ) = PI2ZINV*( FTAN*COSTHE( J+KJ) - FLOG*SINTHE( J+KJ) ) 
GBYNP1(L,J+LJ) = PIZINV*( FLOG*COSTHE(J+KJ) + FTAN*SINTHE( J+KJ)) 
AYNP1(L,J+LJ) = GAYNP1(L,J+LJ)*COSALF(L) -GBYNP1(L,J+LJ)**SINALF(L) 
BYNP1(L,J+LJ) = GAYNP1(L,J+LJ)*SINALF(L)+GBYNP1(L,J+LJ)*COSALF(L) 
CONTINUE 


INFLUENCE COEFFICIENT OF WAKE ELEMENT DUE TO WAKE ELEMENT FROM 
THE OTHER AIRFOIL 


DO 300 L = 1,NAIRFO 

I = NP1*L 

XMID = .5*(X(1) + X(NAIRFO*“NP1+L)) 

YMID = .5*(Y(I) + Y(NAIRFO“NP1+L)) 

KJ = (L-1)*NP1 

MJ = 2*NP] 

NJ = NJ+3 

DXJ = XMID - X(MJ-KJ) 

DXJP = XMID - X(NJ-L) 

DYJ = YMID - Y(MJ-KJ) 

DYJP = YMID - Y(NJ-L) 

FLOG = .5*ALOG((DXJP*DXJP+DYJP*DYJP) /(DXJ*DXJ+DYJ*DYJ) ) 

FTAN = ATAN2(DYJP*DXJ-DXJP**DYJ , DXJP*DXJ+DYJP*DYJ) 

CTIMTJ = COSTHE(1)*COSTHE(MJ-KJ) + SINTHE(1)*SINTHE(MJ-KJ) 
STIMTJ = SINTHE(1)*COSTHE(MJ-KJ) - COSTHE(1)*SINTHE(MJ-KJ) 
GAYNP1(L,NP3) = PI2INV*(FTAN*COSTHE(MJ-KJ) - FLOG*SINTHE(MJ-KJ)) 
GBYNP1(L,NP3) = PI2INV*(FLOG*COSTHE(MJ-KJ) + FTAN*SINTHE(MJ-KJ)) 
AYNP1(L,NP3) = GAYNP1(L,NP3)*COSALF(L)-GBYNP1(L,NP3)*SINALF(L) 
BYNP1(L,NP3) = GAYNP1(L,NP3)*SINALF(L)+GBYNP1(L,NP3)*COSALF(L) 
CONTINUE 
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PIR 1 Be, 


145 
140 
160 


230 
240 
S10, 


INFLUENCE COEFFICIENT ON THE AIRFOIL BY THE WAKE ELEMENT. 


DO 160 L = 1,NAIRFO 
LI=( L-1)*NODTOT 
KI=(L-1)*NP1 

DO 140 I = 1,NODTOT 


XMID = .5*(X(I+KI) + X(I+KI+1)) 

YMID = .5*(Y(I+KI) + Y(I+KI+1)) 

DO 145 MM = 1, NAIRFO 

J = NP1*MM 

DXJ = XMID - X(J) 

DXJP = XMID - X(NAIRFO*NP1+MM) 

DYJ = YMID - Y(J) 

DYJP § = YMID - Y(NAIRFO*NP1+MM) 

FLOG = .5*ALOG((DXJP*DXJP+DYJP*DYJP) /(DXJ*DXJ+DYJ*DYJ) ) 
FTAN = ATAN2(DYJP*DXJ-DXJP*DYJ , DXJP*DXJ+DYJP*DYJ) 
CTIMTJ = COSTHE(I+KI)*COSTHE(J) + SINTHE(1+KI)*SINTHE( J) 
STIMTJ = SINTHE(I+KI)*COSTHE(J) - COSTHE( I+KI)*SINTHE(J) 


AANCI+LI ,2*NODTOT+MM) = PI2ZINV*(FTAN*CTIMTJ + FLOG*STIMTJ) 
BBNCI+LI,2*NODTOT+MM) = PI2INV*(FLOG*CTIMTJ - FITAN*STIMTJ) 
CONTINUE 
CONTINUE 
CONTINUE 


IF GigE) RETR 
MM1 =ehre <1 


INFLUENCE COEFFICIENT ON THE WAKE ELEMENT BY THE CORE VORTICES. 
DO 350 L = 1,NAIRFO 
XMID 0. 5S*CXCNP1*L) + XCNP1*NAIRFO+L) ) 


SCs, On S*CYCNPI*],) 4 Y CN EARN Oto) 
DO 240 MM = 1, NAIRFO 


KN = (MM-1)*MM1 

DO 230 N = 1,MM1 

DX = XMID - XC(MM,N) 
DY = YMID - YC(MM,N) 
DIST2. = DX*DX+DY*"DY 
GCYNP1 = -PI2INV*DX/DIST2 
GCXNP1 = +PI2INV*DY/DIST2 


CYNP1(L,N+KN) = GCYNP1*COSALF(L)+GCXNP1*SINALFCL) 
CXNP1(L,N+KN) = -GCYNP1*SINALF(L)+GCXNP1*COSALF(L) 
CONTINUE 

CONTINUE 

CONTINUE 

IF (NITR. GT.0) RETURN 


INFLUENCE COEFFICIENT ON THE AIRFOIL BY THE CORE VORTICES 


DO 400 L = 1,NAIRFO 


LI =(L-1)*NODTOT 

KI =(L-1)*NP1 

DO 220 I = 1,NODTOT 

XMID = "0. 5*(XCI+KD) xia) 


0. S*CYCI+KI] ee il eee 
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SUMCCN( I+LI ) 0.0 
SUMCCT(C I+LI) 0.0 
DO 280 MM = 1, NAIRFO 


KN = (MM-1)*MM1 

ee 2107 °N = 1,MM1 

Dx = XMID - XC(MM,N) 

DY = YMID - YC(MM,N) 

DIST = SQRT( DX*DX+DY*DY ) 

COSTHN = DX/DIST 

oENEHN = DY/DIST 

CTIMTN = COSTHE(I+KI)*COSTHN + SINTHE( I+KI )*SINTHN 
STIMTN = SINTHECI+KI)*COSTHN - COSTHE( 1+KI)*SINTHN 
CCN =—) *Giirin, Bist 

CCT = -STIMIN/DIST 


SUMCCNCI+LI) = SUMCCNCI+LI) + CCN*CV(MM,N) 
SUMCCT(I+LI) = SUMCCT(CI+LI) + CCT*CV(MM,N) 
210 CONTINUE 
280 CONTINUE 
SUMCCNCI+LI) 
SUMCCT( I+LI) 
220 CONTINUE 
400 CONTINUE 
RETURN 
ERD 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


PI2INV*SUMCCN( I+LI) 
PI2ZINV*SUMCCT( I+LI) 


SUBROUTINE COEF (SINALF,COSALF , OMEGA, UX ,UY,NITR) 


C 
C 
C 
SET COEFFICIENTS OF N EQUS ARISING FROM FLOW C 
TANGENCY CONDITIONS AT MID POINTS OF PANELS C 
SOLVING THE N-SOURCE STRENGTHS IN TERMS OF THE C 
VORTICITY STRENGTH (RESULTING IN 2 RHS) C 
KUTTA CONDITION IS SATISFIED SEPARATELY TO OBTAIN C 
THES VORTICITY STRENGTH C 

THIS SOLUTION METHOD IS DESIRED FOR UNSTEADY FLOW C 
C 

C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


PANEL INDUCED BY UNIT SOURCE DIST ON THE J-TH PANEL 

moljeeso) .... FMRST RHS, COMPONENT OF NORMAL VELO AT THE I-TH 
PANEL WHICH IS DEPENDENT ON VORTICITY STRENGTH OF 
THE FIRS! ALRFOIL 

ACI,NP4) .... SECOND RHS, COMPONENT OF NORMAL VELO AT THE I-TH 
PANEL WHICH IS DEPENDENT ON VORTICITY STRENGTH OF 
THE SECOND AIRFOIL 

Piggies) .... JHIRD RHS, COMPONENT OF NORMAL VELO AT THE I-TH 
PANEL WHICH IS INDEPENDENT OF THE CIRCULATION OF 


C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C iQ). Dee LHS, NORHAL VELOCITY AI THE SE DPOINT «OF THE I-TH 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C BOTH AIRFOILS 

C 


SUBROUTINE COEF (NITR) 

COMMON /BOD/ IFLAG,NLOWER,NUPPER ,NODTOT,X(202),Y(202), 
COSTHE( 201) ,SINTHE( 201) ,SS(2) ,NP1,NP2,NP3,NP4, 
NPS ,XSHIFT, YSHIFT ,NAIRFO,X1( 202) ,Y1( 202), 
COSTHL( 201) , SINTHL( 201) 


t+ + 


ie 


ClaCing 


aqQqona 


C2 Gio) 


ae) 


310 


oA 


ee 
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COMMON /COF/ AC201,211),KEQNS 

COMMON /SING/ Q( 200) ,GAMMA( 2) ,QK( 200) , GAMK( 2) 

COMMON /WAK/ VYW(2) ,VAW(2) ,WAKE( 2) , DT 

COMMON /CORV/ CV(2,200),XC(2,200) ,YC(2, 200) einem yy xt 2) 200), 
CVVY(2,200) ,ACI( 2,200) yYorCza zon, 

COMMON /INF1/ AANC201,201),BBN( 2015201) (AYNETG@22201 ) (By NEw 2920s 


+ SUMAAN( 201,2),SUMBBN( 201, 2) 


COMMON /INF2/ SUMCCN( 201) ,SUMCCT( 201) ,CYNP1( 2,200), 


+ CANEV(C 2,200) 


COMMON /GUST/ UG(200),VG(200) ,XGF ,UGUST , VGUST 

COMMON /MATRIX/ AMAT(201,211) 

COMMON /GEOM/ SINALF(2),COSALF(2) ,OMEGA(2) ,UX(2) ,UY(2),PIVOT, 
XPRM , YPRM 

NEQS § = NODTOT* NAIRFO 

IF (NITR .GT. 0) GO TO 91 


SET/RESET LHS MATRIX A(I,J) FOR EACH NEW TIME STEP 
DO 110 I = 1,2*NODTOT 

DO 110 J = 1,2*NODTOT 

ACI,J) = AANCI,J) 

[PES OiSNE. 0 Gore oa 

FILL IN THE RIGHT HAND SIDE FOR STEADY FLOW 


DO 310 L=1,NAIRFO 


10 = (L-1)*NODTOT 

KI = (L-1)*NP1 

DO 310 I = 1,NODTOT 

XMID = 0.5 * (XI(I+KI) + XI(I+KI+1)) 
YMID = 0.5 * (YICI+KI) + YI(I+KI+1)) 


AC I+LI,NP3) 
AC I+LI ,NP4) 
ACI+LI ,NP5) 
RETURN 


-SUMBBN(CI+LI, 1) 
-SUMBBNCI+LI, 2) 
SINTHLCI+KI)*COSALF(L) - COSTHLCI+KI)*SINALF(L) 


Hoth th ui 


FILL IN THE RIGHT HAND SIDE FOR UNSTEADY FLOW 


DO 210 L=1,NAIRFO 


iol = (L-1)*NODTOT 

KI = (L-1)*NP1 

DO 210 I = 1,NODTOT 

XMID = 0.5 * (XI(I+KI) + XICI+KI41)) 

YMID = 0.5 * (YI(I+KI) + YI(I+KI+1)) 

A(I+LI,NP3) = -SUMBBN(I+LI,1) + BBN(I+LI,NP3)*SS(1)/WAKE( 1) 
A(I+LI,NP4) = -SUMBBN(I+LI,2) + BBN(I+LI,NP4)*SS( 2) /WAKE(2) 
A(I+LI,NPS) = -BBN( I+LI ,NP3)*GAMMA( 1)*SS(1)/WAKE(1)-BBN( I+LI ,NP4) 


+°*GAMMA(2)**SS(2) /WAKE(2) + SINTHL(I+KI)* ((1.+UG(I+LI))*COSALF(L)- 
+VG(I+LI)*SINALF(L)+UX(L)) - COSTHL(I+KI)* ((1.+UG(I+LI) ) 

+*S INALF(L)+VG( I+LI )*COSALF( L)+UY(L) )+OMEGA( L)*( YMID**S INTHL( I+KI) 
+ + XMID*COSTHL( I+KI)) 


ADD CORE VORTEX CONTRIBUTION 


IF (M -EOW 1) GOneeze 
ACI+LI,NP5) = ACI+LI,NPS) - SUMCCNCI+LI) 
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210° SEGNTELNUE 


DO 300 I=1,2*NODTOT 
Be 300 6I=1),2NODLOT+3 
o90 feet l,J) = ACI,J) 


RETURN 
END 


ccccccccccccccccccccccccccccccccccccccccceccccccceccccccccccccccccccccccc 


GAMK 
QK 


Bai)... 
js (21))) re 
Boel) oa. 


AA(I,J) . 


j5)/e 0 a 


C202 CI OOOO IO) a a a ae a 


SUBROUTINE KUTTA (ALPHA, SINALF , COSALF , OMEGA, UX ,UY) 


USING KUTTA CONDITION TO DETERMINE VORTICITY 


GAMMA ... 


gqgqnIaann 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


CIRCULATION FOR THE STEADY FLOW CASE AND ALSO 
CIRCULATION FOR THE PREVIOUS TIME STEP FOR UNSTEADY CASE 
CIRCULATION AT THE CURRENT TIME STEP FOR UNSTEADY CASE 
SOURCE STRENGTH AT THE CURRENT TIME STEP 

QKCI) = B1¢C1)*GAMK( 1)+B2(1)*GAMK(2)+B3(1) 

THAT PART OF THE SOURCE STRNGTH WHICH IS INDUCED BY A 
CIRCULATION OF UNIT STRENGTH FROM AIRFOIL 1 

THAT PART OF THE SOURCE STRNGTH WHICH IS INDUCED BY A 
CIRCULATION OF UNIT STRENGTH FROM AIRFOIL 2 

THAT PART OF THE SOURCE STRNGTH WHICH IS INDEPENDENT 
OF THE CIRCULATION FROM BOTH AIRFOILS 

THAT PART OF THE TANGENTIAL VELOCITY AT THE TRAILING 
EDGE PANELS WHICH IS INDUCED BY A CIRCULATION OF UNIT 
STRENGTH BY AIRFOIL J 

THAT PART OF THE TANGENTIAL VELOCITY AT THE TRAILING 
EDGE PANELS WHICH IS INDEPENDENT OF THE CIRCULATION 


SUBROUTINE KUTTACNITR,PVTAG ) 
COMMON /BOD/ IFLAG,NLOWER,NUPPER,NODTOT, X( 202) ,Y¥( 202), 


a 
at 
Be 


COSTHECZ0] ), SINTHE( 201) ~SS(2) ,NP1,NP2,NP3,NP4, 
NMeo{ASHIFI ,YSHIFT,NAIRFU, x10 202) ,YI( 202), 
COSTHL( 201) ,SINTHL( 201} 


COMMON /COF/ AC 201,211) ,KEQNS 

COMMON /SING/ Q( 200) ,GAMMA( 2) ,QK( 200) ,GAMK( 2) 

COMMON /WAK/ VYW(2),VXW(2) ,WAKE(2) , DT 

COMMON /CORV/ CV(2,200) ,XC(2,200) ,YC( 2,200) ,M,TD,CVVX(2,200), 


+ 


Cy xC2. 200 )exemee, 200) erCl(2, 200) 


COMMON /INF1/ AAN(201,201),BBN(201,201) ,AYNP1(2,201),BYNP1(2,201), 


a 


SUMAAN( 201,2) ,SUMBBN( 201, 2) 


Gonmen 7INF2/ SUMCCN(201),SUMGET( 201) ,CYNP1(2,200), 


ne 


CXNP1( 2,200) 


COMMON /GUST/ UG(200) ,VG(200) ,XGF , UGUST, VGUST 

COMMON /MATRIX/ AMAT( 201,211) 

COMMON /BMAT/ B1,B2,B3 

COMMON /GEOM/ SINALF(2),COSALF(2),OMEGA(2),UX(2),UY(2),PIVOT, 


a 


XPRM, YPRM 


COMMON /GIES/ NGIES 

DIMENSION B1(200),B2( 200) ,B3(200) ,AA(4,2),BB(5) ,EE(2),FF(2), 
+GG( 2) ,RADI(2),AAA( 2) , BBB( 2) ,CCC( 2) , DDD(2) ,EEE(2) ,FFF(2), 
+DG(2),E(2),GAMA(2) ,DELG(2) , DGAMK( 2,3) ,COEFL(2,3) 


C RETRIEVE SOLUTION FROM A-MATRIX FOR THE SOURCE STRENGTH IN TERMS 


Oly 


Ci CGC 


2 © qa 


qaqa 


50 


419 
430 


425 


450 


460 


470 


THE TWO GAMK CONTRIBUTION AND THE CONSTANT COMPONENT 


COUNT =0 
DO 50 I = 1,NODTOT*2 
B1(1) = ACI,NP3) 
B2(1) = ACI,NP4) 
B3(I) = ACI,NPS) 


IF (M.GT.0) GOTO 400 
STEADY KUTTA CONDITION 
COMPUTE TANGENTIAL VELOCITIES 


DO 425 L = 1,NAIRFO 

LI = (L-1)*NODTOT 

KI = (ib eNedl 

KK = (L-1)*2 

DO 430 K=1,2 

IF (K -EQ291)31 =e 

IF (K .EQ 2) 1 = NODpTOr 

AMID = 0.5 * (A0GiRD) Glee Ie 1) 
Yo = 0.593 Cok) Seer) 
AAC K+KK,1) = SUMAAN(CI+LI,1) 

AACK+KK,2) = SUMAANCI+LI,2) 

BBC K+KK) = COSALF(L)*COSTHLCI+KI) + SINALFC(L)*SINTHL(I+KI) 


DO 419 J = 1,2*NODTOT 


AAC Kt+KK, 1) = AACK+KK,1) - BBNCI+LI,J)*B1(J) 
AACK+KK, 2) = AACK+KK,2) - BBNCI+LI,J)*B2(J) 
BB( K+KK) = BBC K+KK) SEN Itbt Bs) 
CONTINUE 

CONTINUE 

CONTINUE 


SET UP KUTTA CONDITIONS 


DO 450 I=1,2 

LI = (I-1)*NAIRFO+1 
DGAMK(I,1) = AA(LI,1) + AA(LI+1,1) 
DGAMK(1,2) = BACLI.2) + SACnIe Ee) 
DGAMK(I,3) = -(BB(LI) + BB(LI+1)) 


SOLVE KUTTA CONDITION BY GAUSS 


R=DGAMK( 2,1) /DGAMK(1,1) 

DO 460 K = 2,3 

DGAMK(2,K) = DGAMK(2,K)-R*DGAMK(1,K) 

GAMMA(2) = DGAMK(2,3)/DGAMK(2,2) 

GAMMA(1) = (DGAMK(1,3)-DGAMK( 1, 2)*GAMMA( 2) ) /DGAMK(1, 1) 


CALCULATE SOURCE STRENGTH 


DO 470 L = 1,NAIRFO 


LI =(L-1)*NODTOT 
DO 470 I = 1,NODTOT 
Q( I+LI) = GAMMA(1)*B1CI+LI) + GAMMA(2)*B2(I+LI) + B3(CI+LI1) 


RETURN 
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CG Care) 


Cr <2 62 


CarG2-03 


UNSTEADY KUTTA CONDITION 
FIND VI AT TRAILING EDGE PANELS 


INITIAL GUESS FOR GAMK FROM PREVIOUS TIME STEP 
oO Mme CNITR .EQ. 0) THEN 


GAMK( 1) = GAMMA(1) 

GAMK(2) = GAMMA( 2) 

ENDIF 

DO-125 L = 1,NAIRFO 

LI = (L-1)*NODTOT 

KI = (L-1)*NP1 

KK = (L-1)*2 

mise” K = 1,2 

Mm (K.EQ. 1) 1 = 1 

Peek EQ. 2) 1 = NODITOT 

XMID =0.5 *“S@ atk jatexICI+KI+1)) 
TID = OSes (yi) + YC I+KI+1)) 


VTANG = ((1.4+UG(I+LI))*COSALF(L)-VG(I+LI)*SINALF(L)+UX(L)) 
+**COSTHL( I+KI)+ (( 1. +UG(I+LI) )*SINALF( L)+VG( I+LI )*COSALF( L)+UY(L)) 
+*SINTHL(I+KI) + OMEGA( L)*( YMID*COSTHL( I+KI) 
+- XMID*SINTHL( I+KI)) 

AAC(K+KK,1) = - AAN(I+LI,NP3)*SS(1)/WAKE(1) 
AA(K+KK,2) = - AANCI+LI,NP4)*SS(2)/WAKE(2) 
BB(K+KK) = VTANG + AAN(I+LI,NP3)*SS(1)*GAMMA(1)/WAKE(1) + 
+ AAN( I+LI ,NP4)**SS( 2)*GAMMA( 2) /WAKE( 2) 
DO 119 J = 1,NODTOT 
AA(K+KK,1) = AAC(K+KK,1) + AAN(I+LI,J) - BBN(I+LI,J)*B1(J) 
AA(K+KK,2) = AA(K+KK,2) + AAN(I+LI,J+NODTOT) - BBN(I+LI,J)*B2(J) 
BB(K+KK) = BB(K+KK) - BBN(I+LI,J)*B3(J) 
119 CONTINUE 
DO 120 JJ = NODTOT+1,2*NODTOT 
AA(K+KK,1) = AAC(K+KK,1) - BBNCI+LI,JJ)*B1(JJ) 
AA(K+KK,2) = AAC(K+KK,2) - BBN(I+LI,JJ)*B2(JJ) 
BB(K+KK) = BB(K+KK) - BBN(I+LI,JJ)*B3(JJ) 
120 CONTINUE 


ADD CORE VORTEX CONTRIBUTION 


im (i, Ge. 1) GOTO 100 

BBC K+KK) = BB(K+KK) + SUMCCTCI+LI) 
100 CONTINUE 
130 CONTINUE 
125. CONTINUE 

ie CNGIES .EQ. 1) GOTO 145 


SATISFYING KUTTA CONDITION -- SOLVE FOR VORTEX STRENGTH 


be 135.) = 1,2 


LI = (I-1)*NAIRFO+1 

AAA(I) = AACLI,1)**2-AACLI+1,1)**2 

BBB(I) = AACLI,2)**2-AACLI+1,2)**2 

Geet = 2 Chr L1,1)*BBC LI ) -AA( L1+log) )*BBC LI+1) -(2-1)*S5(1)/DT) 
DDD(I) = 2*(AACLI,2)*BBCLI)-AAC LI+1 ,2)*BBC(LI+1)-(1-1)*SS(2)/DT) 
EEE(I) = 2*(AACLI,1)*AACLI,2)-AACLI+1,1)*AACLI+1, 2) ) 


rrr (1) BB( LI )***2-BB( LI+1)**2+2*SS( I )*GAMMA( I) /DT 
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135 CONTINUE 


60 DO 140 I = 1,2 
DGAMK(1I,1) = 2*AAACT)*GAMK(1)+EEE( I )*GAMK( 2)+CCC( I) 
DGAMK(1I,2) = 2*BBB( I)*GAMK( 2)+EEE(1)*GAMK(1)+DDD( I) 


140 DGAMK(I,3) = -CAAACI)*GAMK( 1)**2+CCC(1I)*GAMK(1)+BBB(1)* 
+ GAMK(2)%**2+DDD( I )*GAMK( 2)+EEE( I)*GAMK( 1)*GAMK( 2)+FFF(1)) 
C140 WRITE (6,*)'DGAMK' ,I,DGAMK(1I,1),DGAMK(I,2),DGAMK(I,3) 
DO 144 I=1,2 
DO 144 KKK=1,3 
144 COEFL(I,KKK) = DGAMK(I, KKK) 


GAUSSIAN ELIMINATION 


ca Gir@ 


R=DGAMK( 2,1) /DGAMK(1, 1) 
DO 200 K = 2,3 
200 DGAMK(2,K) = DGAMK(2,K)-R*DGAMK(1,K) 
DG(2) = DGAMK(2,3)/DGAMK( 2,2) 
DG(1) = (DGAMK(1,3)-DGAMK(1,2)*DG(2)) /DGAMK(1,1) 
GAMK(1) = GAMK(1)+DG(1) 
GAMK(2) = GAMK(2)+DG(2) 
IF ((ABS(DG(1)/GAMK(1)) .LT..0001) . AND. 
+ (ABS(DG(2)/GAMK(2)).LT..0001)) GO TO 300 
COUNT = COUNT + 1 
IF (COUNT .EQ. 50) GO TO 310 
GO TO 60 
310 WRITE(6,*)'KUTTA CONDITION NOT SATISFIED -- COUNT EXCEEDED 50' 
RETURN 
C 
C SET UP KUTTA CONDITIONS FOR GIESING CONDITION 
C 
145 DO 451 I=1,2 
il = (I-1)*NAIRFO+1 
DGAMK(I,1) = AA(LI,1) + AA(LI+1,1) 
DGAMK(I,2) = AA(LI,2) + AA(LI+1,2) 
451 DGAMK(I,3) = -(BB(LI) + BB(LI+1)) 
fc 
C SOLVE KUTTA CONDITION BY GAUSS 
C 


R=DGAMK( 2,1) /DGAMK(1, 1) 
DO 461 K = 2,3 
461 DGAMK(2,K) = DGAMK(2,K)-R*DGAMK(1,K) 
GAMK(2) = DGAMK(2,3)/DGAMK(2, 2) 
GAMK(1) = (DGAMK(1,3)-DGAMK( 1, 2)*GAMK( 2) )/DGAMK(1,1) 


CALCULATE SOURCE STRENGTH 


aaa 


300 DO 160 L = 1,NAIRFO 
LI =(L-1)*NODTOT 
DO 160 I = 1,NODTOT 
160 QKCI+LI) = GAMK(1)*B1(I+LI) + GAMK(2)*B2(I+LI) + B3(I+LI) 


CALCULATE TANGENTIAL VELOCITY AT THE TRAILING EDGE BY BACK 
SUBSTITUTION. THE PRODUCT OF THE TRAILING EDGE VELOCITY 
SHOULD BE NEGATIVE . THIS IS CHECKED IN THE MAIN PROGRAM. 


C2 Canc) Ca C2 


VTAG1 = AA(1,1)*GAMK(1)+AA(1,2)*GAMK(2)+BB(1) 


(AGA IS OL 


VTAG5O AA(2,1)*GAMK( 1)+AA( 2, 2)*GAMK( 2)+BB( 2) 
VTAG51 AA(3,1)*GAMK( 1)+AA( 3, 2)*GAMK(2)+BB(3) 
VTA100 = AA(4,1)*GAMK(1)+AA(4,2)*GAMK(2)+BB(4) 
PVTAG = VTAG1*VTAGS5O 

FORMAT( 1X, 'VTAG' ,4F10. 6) 

RETURN 

END 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


C2 C2 C2Caec? C2: Ga: Ca Garg C2 C3 
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+++ 


SUBROUTINE TEWAK (SINALF,COSALF) 
COMPUTE WAKE ELEMENT AT THE TRAILING EDGE 


CQ) 62 Co) C2 


WOWKCL) .... Y-VELO COMP AT THE MID-POINT OF THE WAKE ELEMENT 
WITH RESPECT TO THE LOCAL FROZEN FRAME OF REFERENCE 

Verne) .....X-VELO COMP AT THE MID-POINT OF THE WAKE ELEMENT 
WITH RESPECT TO THE LOCAL FROZEN FRAME OF REFERENCE 


SUBROUTINE TEWAK 


COMMON /BOD/ IFLAG,NLOWER ,NUPPER,NODTOT,X( 202) ,Y(202), 
COSTHE( 201) ,SINTHE( 201) ,SS(2) ,NP1,NP2,NP3,NP4, 
NP5 ,XSHIFT, YSHIFT ,NAIRFO,X1(202) ,YI(202), 
COSTHL( 201) , SINTHL( 201) 
COMMON /COF/ A(201,211),KEQNS 
COMMON /SING/ Q( 200) ,GAMMA( 2) ,QK( 200) ,GAMK( 2) 
COMMON /WAK/ VYW(2),VXW(2) ,WAKE(2) ,DT 
COMMON /WAK2/ VYWK(2),VXWK(2) 
COMMON /CORV/ CV(2,200) ,XC(2,200),YC(2,200) ,M,TD,CVVX(2,200), 
OvVy02.200) ,.XCI( 2.200) . YCI(2, 200) 
COMMON /INF1/ AAN(201,201),BBN(201,201),AYNP1(2,201),BYNP1(2,201), 
SUMAAN( 201, 2) , SUMBBN( 201, 2) 
COMMON /INF2/ SUMCCN( 201) ,SUMCCT(201),CYNP1(2,200), 
CXNP1( 2,200) 
COMMON /GEOM/ SINALF(2),COSALF(2) ,OMEGA(2) ,UX(2) ,UY(2) ,PIVOT, 
XPRM, YPRM 
COMMON /GUST/ UG( 200) ,VG(200),XGF ,UGUST, VGUST 
DIMENSION CONT1X(2) ,CONT1Y(2) ,CONT2X(2) ,CONT2Y(2), 


+ CONT3X(2) ,CONT3Y(2) ,CONT4X(2) ,CONT4Y(2) ,CONT5X( 2) , CONT5Y( 2) 


DO 20 I 
NN 

XMID 
YMID 
UGW 
VGW 0.0 

XG XMID 

Te (3G .GI.. AGF) GO TO 10 

UGW UGUST 

VGW VGUST 

VYwK(1) (1. +UGW)*SINALF(1)+VGW*COSALF( TI) 
VXWK(T) (1. +UGW)*COSALF(I) -VGW*SINALF( TI) 
FORMAT(1X,'SINALF :' ,4F10. 5) 


1,NAIRFO 

(I-1)*NP1 

0.5 * (X(NP1+NN) + X(NP24I-1)) 
0.5 * (Y(NP1+NN) + Y(NP24I-1)) 
0.0 


CONT1Y(I) = (1.+UGW)*SINALF( I )+VGW*COSALF(I) 
CONT1X(1) = (1.+UGW)*COSALF(I ) -VGW*S INALF(1) 
Commer) = 0 


12] 


C262. C2 


P20 
110 


© 2 © 


is 


C 
C 
C 


130 
150 
140 
20 
Pb Jag 
2225 


CONT2X( 1) 
CONT3Y( 1) 


0 
0 
CONT3X(T) 0 


CONTRIBUTION FROM THE SOURCE AND WAKE DIST FROM THE AIRFOIL 


DO 110 L = 1,NAIRFO 

KJ = (L-1)*NODTOT 

DO 120 J = 1,NODTOT 

VYWK(I) = VYWK(I) + AYNP1(1,J+KJ)*QK(J+KJ) + BYNP1(1I,J+KJ)*GAMK(L) 
VXWK(I) = VXWK(I) - BYNP1(1,J+KJ)*QK(J+KJ) + AYNP1(1,J+KJ)*GAMK(L) 
CONT2Y(1) = CONT2Y(1)+ AYNP1(1,J+KJ)*QK(J+KJ) 

CONT2X(I) = CONT2X(I)- BYNP1(1,J+KJ)*QK( J+KJ) 

CONT3Y(I) = CONT3Y(I)+ BYNP1(1,J+KJ)*GAMK(L) 

CONT3X(1) = CONT3X(I)+ AYNP1(1I,J+KJ)*GAMK(L) 

CONTINUE 


ADD WAKE ELEMENT CONTRIBUTION OF THE OTHER AIRFOIL 


J = 3-I 
VYWK(I) = VYWKCI)+ BYNP1(1,NP3)*( GAMMA( J) -GAMK(J)) 


+ * S5(J) /WAKE(CJ) 


VXWK( I) = VAWK(1I)+ AYNP1(1,NP3)*( GAMMA( J) -GAMK(J)) 


+ * SS(J)/WAKE(J) 


CONT4Y( I) = BYNP1(1,NP3)*(GAMMA( J) -GAMK( J) ) 


+ * SS(J)/WAKE( J) 


CONT4X(1) = AYNP1(1,NP3)%*(GAMMA(J)-GAMK(J)) 


+ * SS(J)/WAKE( J) 


ADD CORE VORTEX CONTRIBUTION 


CONTS5Y(I) = 0 
CONTSX(I) = 0O 

IF (M .EQ. 1) GO TO 140 
MM1 =e 

DO 150 L = 1,NAIRFO 

KN =(L-1)*MM1 


DO 130 N = 1,MM1 

VYWK(I) = VYWK(I) + CYNP1(1,N+KN)*CV(L,N) 
VXWK(1) = VXWK(I) + CXNP1(I,N+KN)*CVCL,N) 
CONT5Y(I) = CYNP1(1I,N+KN)*CV(L,N) 
CONT5X(I) CXNP1( I ,N+KN)*CV(L,N) 
CONTINUE 

CONTINUE 

CONTINUE 

FORMAT (1X, 'VXWK :',6F10.5) 

FORMAT (1X, 'VYWK :',6F10. 5) 

RETURN 

END 
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SUBROUTINE PRESS (SINALF,COSALF ,OMEGA,UX,UY) 


COMPUTE UNSTEADY FLOW PRESSURE DISTRIBUTION 
AND VELOCITY POTENTIAL AT MID-POINTS OF PANELS 


Cart C26) Gag 


eZ 


cccccccccccccccccccccecccccccccccccccccccccccccccccccccccccccccccccccccc 


QqOGe AON am oI 0) Orr a Oo rea Oe OO a a) ee es ee 


WIAMNG ..... TANGENTIAL VELO COMP OF A FIXED POINT OF THE MOVING 
LOCAL FRAME OF REFERENCE. 

PHOKGI)...<. DISTURBANCE VELO POTENTIAL AT THE MID POINT OF THE 
I-TH PANE AT THE CURRENT TIME STEP FOR THE UNSTEADY 
FLOW CASE. 

ree) ..... DISTURBANCE VELO POTENTIAL AT THE MID POINT OF THE 
PS iieraNEt POR ane PREVIOUS TIME STEP 

PHILE(L).... DIFFERENCE OF THE POTENTIALS OF THE LEADING EDGE TO 
THE LOWER TRAILING EDGE FOR THE RESPECTIVE AIRFOIL 

PINK(L) .... DIFFERENCE OF THE POTENTIALS AT A POINT 100 CHORD 
LENGTH UPSTREAM OF THE LE FOR THE RESPECTIVE AIRFOIL 

SUMC(L) .... GAMMA ASSOCIATED WITH THE INTEGRATION OF THE DISTURB- 


ANCE VELOCITY AROUND THE WHOLE AIRFOIL 


THE FOLLOWING INFLUENCE COEF ARE COMPUTED FOR A FINER GRID ON THE 
AIRFOIL SO AS TO OBTAIN A MORE ACCURATE VELO POTENTIAL AT THE TE 
THE INFLUENCE COEF ON THE I-TH PANEL FROM THE J-TH PANEL OF THE 
AIRFOIL REMAINS THE SAME FOR ALL TIME STEP 


AANP1(1I,J,K).. NORMAL VELOCITY INDUCED AT THE I-TH PANEL SUB NODE K 


OF AIRFOIL 1 DUE TO UNIT STRENGTH DIST SOURCE STRENG- 
THON THE J-TH PANEL OF AIRFOIL 1 


BBMP1(1,J,K).. NORMAL VELOCITY INDUCED AT THE I-TH PANEL SUB NODE K 


OF AIRFOIL 1 DUE TO UNIT STRENGTH DIST VORTICITY STR- 
ENGUHO ON THE J-fH PANEL OF VAIRFOIL 1 


AANP2(1,J,K).. NORMAL VELOCITY INDUCED AT THE I-TH PANEL SUB NODE K 


OF AIRFOIL 2 DUE TO UNIT STRENGTH DIST SOURCE STRENG- 
PaO tne J-ii PANEL OF AIKRFORD 2 


BBNP2(I,J,K).. NORMAL VELOCITY INDUCED AT THE I-TH PANEL SUB NODE K 


+ 
+ 
cL 


Ee 
4 


+ 


OF VAUREOIG 2 DUE TO UNIT SsaheNGih, DIST YORTICITY SIR- 
ExGiHeOh THE J-iH PANEL OF AIRFOIL 2 


SUBROUTINE PRESS 


COMMON 


COMMON 
COMMON 
COMMON 
COMMON 
COMMON 


COMMON 
COMMON 
COMMON 


COMMON 
COMMON 


/BOD/ IFLAG,NLOWER ,NUPPER ,NODTOT, X(202) ,Y(202), 
COSTHE( 201) , SINTHE( 201) ,SS(2),NP1,NP2,NP3,NP4, 
NPS ,XSHIFT , YSHIFT ,NAIRFO,X1I(202) ,YI(202), 
COSTHL( 201) , SINTHL( 201) 
/CPD/ CP( 200) ,SCL,T,SCM,SGAM 
/NUM/ PI,PI2INV 
/SING/ Q( 200) ,GAMMA( 2) ,QK( 200) ,GAMK(2) 
/WAK/ VYW( 2) ,VXW(2) ,WAKE(2) ,DT 
/CORV/ CV(2,200) ,XC( 2,200) ,YC(2,200) ,M,TD,CVVX(2,200), 
CVVY( 2,200) ,XC1(2, 200) , YCI(2,200) 
/INF1/ AAN( 201,201) ,BBN( 201,201) ,AYNP1(2,201) ,BYNP1(2,201), 
SUMAAN( 201,2) , SUMBBN( 201, 2) 
/INF2/ SUMCCN( 201) ,SUMCCT( 201) ,CYNP1( 2,200), 
CXNP1( 2,200) 
/INFL3/ AANP(101,101,6),BBNP(101,101,6) 
/POT/ PHI( 200), PHIK( 200) 
/GUST/ UG( 200) ,VG( 200) ,XGF,UGUST, VGUST 
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qOoaa an 


COMMON /EXTV/ UE(200),VN(200) 
COMMON /BMAT/ B1,B2,B3 
COMMON /GEOM/ SINALF(2),COSALF(2),OMEGA(2) ,UX(2),UY(2),PIVOT, 
“ XPRM , YPRM 
DIMENSION PHITEL(2),PHITEU(2),WGHT(5),PLOC(5),SUMC(2), 
+AANP1(50,50,9),AANP2(50,50,9), BBNP1(50,50,9), BBNP2(50,50,9), 
+COSTHP( 102,6) , SINTHP( 102,6),AANP4(2),PHIL( 102) ,UGU(100,6), 
+VGU(100,6) ,PDUM(2) ,DP(2) , BBNP4(2) ,B1(200) ,B2( 200) ,B3( 200) , PHILE(2) 
3333 DIMENSION CON2(2),VIGT(2,50) ,WAKCON( 2,50), VOTCON( 2,50), 
rm DUMCON( 2,50) 
7733 DIMENSION CP1(200) 
4444 REAL *8 PINK(2),VELX 
DATA WGHT/. 11846344,. 23931434,. 28444444, 


4 .23931434,.11846344/ 
DATA PLOC/. 04691008,. 23076535,. 50000000, 
+ . /6923466,. 95300679/ 


WRITE (6, 1000) 
IF (NM .GT. 1) GO TO 510 


COMPUTE THE INFLUENCE COEFFICIENT FOR A FINER GRID 
ON THE AIRFOIL BY THE SAME AIRFOIL ... 
( ONLY COMPUTED ONCE ) 
DO 200 L = 1,NAIRFO 
LI = (L-1)*NODTOT 
KI = (L-1)*NP1 
DO 200 I = 1,NODTOT 
DO 200 K =1,5 
DX = PLOC(K)*(X( I+KI+1)-X(1+KI)) 
DY = PLOC(K)*(Y(I+KI+1)-Y(I+KI)) 
XMID = X(I+KI) + DX 
YMID = Y(I+KI) + DY 
DO 200 J = 1,NODTOT 
FLOG = 0.0 
FTAN = PI 
IF (J+LI .EQ. I+LI) GO TO 100 
DXJ = XMID - X(J4KI) 
DXJP = XMID - X(J+KI+41) 
DYJ = YMID - Y(J+KI) 
DYJP § = YMID - Y(J+KI41) 
FLOG = .5*ALOG((DXJP**DXJP+DYJP*DYJP ) /(DXJ*DXJ+DYJ*DYJ) ) 
FTAN = ATAN2(DYJP*DXJ-DXJP*“DYJ , DXJP*DXJ+DYJP*DYJ) 
100 CTIMTJ = COSTHE(I+KI)*COSTHE(J+KI) + 
s. S INTHE( I+KI )*SINTHE(J+KI) 
STIMTJ = SINTHE(I+K1)*COSTHE(J+KI) - 
4. COSTHE( 1+K1)*SINTHE( J+KI ) 
IF G.EQ. 2) Gonos 
AANP1(I,J,K) = PIZINV*(FTAN“CTIMTJ + FLOG*STIMTJ) 
BBNP1(I,J,K) = PIZINV*(FLOG*CTIMTJ - FTAN*STIMTJ) 
GOTO 200 
120 AANP2(I,J,K) = PIZINV*(FTAN*CTIMTJ + FLOG*STIMTJ) 
BBNP2(1,J,K) = PI2ZINV*(FLOG*CTIMTJ - FTAN*STIMTJ) 


200 CONTINUE 
510 CONTINUE 
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FIND TANGENTIAL VELOCITY AT THE MIDPOINT OF THE I-TH PANEL 


DO 600 L = 1,NAIRFO 
LI=(L-1)*NODTOT 
KI=(L-1)*NP1 


SUMC(L) =0.0 

DO 600 I = 1,NODTOT 

CONTR1 = oC. U 

CONTR2 = 0.0 

WAKCON(L,I) = 0.0 

VOTCON(L,I) = 0.0 

DUMCON(L,I) = 0.0 

XMID = 07S te AICIERD) + XICI+KI41)) 
YMID = 0.5 % (YICI+h8) + YI@HKI+41)) 
DX = (X1GR4KI41) see XTC 1#KI)) 

DY = (yU@ie tl) = Yt 14+KT)) 

ist = SQRT(DX*DX+DY*DY ) 


ACCOUNT FOR THE FREESTREAM AND THE MOTION OF THE AIRFOIL 


Vox 
4 


WSY = 


(1. +UG(I+LI) )*COSALF(L)-VG(I+LI)*SINALF(L) + OMEGA(L)* 
¥MID + UX(L) 
(1. +UGC I+L1) )*SINALF(L)+VG( I+LI)*COSALF(L) 


+ - “OMEGA L)“ALD + UYCL) 
VS = VSH*VSX + WSY>VSY 


VNORM = 
VIANG 


-VSR=SINTHL( I+KI )#¥VSY*COSTHL(1I+KI) 


= (((1. +UG( I+LI) )*COSALF(L)-VGC I+LI)*SINALF(L)+UX(L)) 


eee UCIT aR + (C1 +UG( ISLI) )*SINALF(L)+VGCI+LI1) *COSALF( L)+UY(L)) 
= “oUNGHE( (+R1)+ OMEGACL)*(YMID*COSTHL(I+KI) 


+  - XMID*SINTHL(I+KI))) 


VOPREE = 
VACT = 


INTRODUCE 
VELO ONLY 


DO 260 
DX 

Ue 
ial 
YIN’ 
VDUM 


nA 


niu uw u 


INFLUENCE 


DO 245 MM 
=) 

DXJ 
DXJP 
bx J 

bi PP 
FLOG 
FTAN 
CT J 
4. 

od iis 5 Ji] 


VTANG 
VTANG 


oaiginbee Rt DS FOR THE PURPOSE OF “THE VELO PODSNITIAL. 
CALCULATED AT THE MIDPT OF THE PANEL WHERE K = 


= 1,5 
PEOGCe)-(XC1+KiI+1)-X( ER1)) 
PLOC(K)*CY¥CI+KI+1)-YCI+KI)) 

X(I+KI) + DX 

YC I+KI) + DY 

0.0 


COEF ON I-TH PANEL DUE TO THE WAKE ELEMENTS 


= 1, NAIRFO 

NP 1*MM 

Kiw - XCJ) 

XINT - XCNAIRFO*NP1+MM) 

Wis - Y(J) 

YIN - YCNAIRPO*NP1+11) 

. 5*ALOG( (DXJP*DXJP+DYJP*DYJP) /(DXJ*DXJ+DYJ*DYJ) ) 
ATAN2 (DY JP*DXJ-DXJP*DYJ , DXJP*DXJI+DYJIP*DYJ) 
COSTHECI+KI)*COSTHE(J) + 

S.INTHE ( I+K1)*S2NTHE( J) 

Sitene( I+KI)*COSTHECJ) - 


245 


Cr Clic 


oC). C2 <2 C2 Go 
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270 


+ COSTHE( I+KI)*SINTHE( J) 
AANP4(MM) = PI2ZINV*(FTAN*CTIMTJ + FLOG*STIMTJ) 
BBNP4(MM) = PI2ZINV*(FLOG*CTIMTJ - FTAN*STIMTJ) 
CONTR1 = SS(1)*(GAMMA( 1) -GAMK( 1) )*AANP4( 1) 
+/WAKE(1)+SS(2)*(GAMMA( 2) -GAMK( 2) )*AANP4( 2) /WAKE( 2) 


CONTRIBUTION TO VELO COMPONENT BY WAKE ELEMENT. 


IF (Kk . EQjgs) ene. 

VACT = VACT + SS(1)*(GAMMA( 1)-GAMK(1))*AANCI+LI, 
+NP3) /WAKE(1)+SS(2)**( GAMMA( 2 ) -GAMK( 2) )*AANCI+LI,NP4) /WAKE( 2) 
VNORM = VNORM + SS(1)**(GAMMA( 1) -GAMK( 1) )**BBNP4( 1) 
+/WAKE(1)+SS(2)*( GAMMA( 2) -GAMK( 2) )*BBNP4( 2) /WAKE( 2) 

ENDIF 


EFFECTS ON AIRFOIL BY THE SAME AIRFOIL 


INTEGRATION AROUND FIRST AIRFOIL 
CONTRIBUTION TO VELO COMP BY AIRFOIL 1 WHEN K = 3 


DO 300° J°= 1 Nerier 

IF (L .EQ.2 2 ec0te rg 

VDUM = VDUM-BBNP1(1,J,K)*QK(J)+AANP1(1,J,K)*GAMK( 1) 
IF (K .EQ. 3) 

VACT = VACT-BBNCI+LI,J)*QK(J)+ 
+AANCI+LI , J)*GAMK( 1) 

VNORM = VNORM+( AANP1(1,J,K)*QK( J) )+ 
+( BBNPIC1,J,K)*GAMNK(1)) 

ENDIF 

GOTO 300 


INTEGRATION AROUND SECOND AIRFOIL 
CONTRIBUTION TO VELO COMP BY AIRFOIL 2 WHEN K = 3 


VDUM = VDUM-BBNP2(1I,J,K)*QK(J+NODTOT)+AANP2(I,J,K)**GAMK(2) 
IF (K .EQ. 3) THEN 

VACT = VACT-BBN( I+LI , JtNODTOT)*QK( JtNODTOT)+ 
+ AAN( I+LI , J+NCDTOT )*GAMK(2) 

VNORM = VNORM+( AANP2(1,J,K)*QK( J+NODTOT) )+ 
+ (BBNP2(I,J,K)*GAMK( 2) ) 

ENDIF 


300 CONTINUE 
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EFFECTS ON AIRFOIL BY THE OTHER AIRFOIL 


INTEGRATION AROUND BOTH AIRFOIL 
CONTRIBUTION TO VELO COMP BY BOTH AIRFOIL WHEN K = 3 


DO 350 J = NODTOT+2,2*NODTOT+1 

DXJ = XINT - X(J-KI) 

DXJP = XINT - X(J-KI+1) 

DYJ = YINT - Y(J-KI) 

DYJP = YINT - Y(J-KI+1) 

FLOG = .5*ALOG((DXJP*“DXJP+DYJP*DYJP) /( DXJ*DXJ+DYJ*DYJ) ) 


oe kekere) 


350 


400 


260 


11S) 


600 


FTAN ATANGCDY JR Wx Deas J, DX JP*DXI+DYIP*DY J ) 


CTIMTJ COSTHE( I+KI )*COSTHE(J-KI) + 
+ SINTHE( I+KI)*SINTHE( J-KI) 

Sry == SINTHEQCIA! )eeslhEt=K1) - 
+ COSTHE( I+KI )*SINTHE( J-KI) 

AANP3 Pret Chan Geli + FLOG*STIMTJ) 


BBNP3 = PI2ZINV*(FLOG*CTIMTJ - FTAN*STIMTJ) 
ee Ko.) THEN 
VACT = VACT-BBN( I+LI , J-LI-1)*QK( J-LI-1)+ 


+ AAN( I+L1,J-LL-1)*GAMK(3-L) 

VNORM = VNORM+( AANP3*QK( J-LI-1))+ 
+ ( BBNP3*GAMK( 3-L) ) 

ENDIF 

VDUM = VDUM - BBNP3*QK(J-LI-1) + AANP3*GAMK(3-L) 


CONTRIBUTION BY CORE VORTICES 
1. TO VELO POTENTIAL 
2. TO VELO COMPONENT ONLY WHEN K = 3 


Peeen- EQ. 1) GOTO 150 


MM1 = eve? 

SUMCN = 0.0 

SUMCT = 0.0 

DO 400 MM = 1, NAIRFO 

KN = (MM-1)*MM1 

Poes00 Ne= 1 Ml 

DXC = Adel = ACCOM SN) 

DYC = YINT - YCCMM,N) 

DISTC = SQRT(DXC*DXC+DYC*DYC) 

COSTHN = DXC/DISTC 

SINTHN = DYC/DISTC 

Clive = COSTHE(I+KI)*COSTHN + SINTHE(I+KI)*SINTHN 
STIMTN = SINTHECI+KI)*COSTHN - COSTHE( I+KI)*SINTHN 
CCN = -CTIMIN/DISTC 

CCT = -STIMIN/DISTC 

SUMCN = SUMCN + CCN*CV(MM,N) 

SUMCT = SUMCT + CCT*CV(MM,N) 


eet kK . £Q5 3) THEN 
VACT = VACT+SUMCT*PI2INV 
VNORM = VNORM+( PI2ZINV*SUMCN) 


ENDIF 

CONTR2 = PI2INV*SUMCT 

CONTINUE 

VTIANG VTANG + DBLE( (CONTR1+CONTR2+VDUM )*WGHT(K) ) 


SUMC(L) = SUMC(L)+VDUM*DIST*WGHT(K) 
WAKCON(L,I) = WAKCON(L,1I)+CONTR1 


VOTCON(L,I) = VOTCON(L,1I)+CONTR2 
DUMCON(L,I) = DUMCON(L,1I)+VDUM 
CONTINUE 

PHIK(I+LI) = ((VTANG)-VTFREE)*DIST 
CP(I+LI) = VS - (VACT*VACT) 
CP1(I+LI) = CP(I+LI) 

UE(I+LI) = VACT 

VN(I+LI) = VNORM 

CONTINUE 


6668 FORMIATCZI5,7E15.7) 
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3335 FORMAT (215,7F10. 6) 
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COMPUTE DISTURBANCE POTENTIAL AT THE LEADING EDGE BY LINE 
INTEGRAL OF THE VELOCITY FIELD 
FROM UPSTREAM (AT INFINITY) TO THE LEADING EDGE 


TRY =} 

PINK(1) = 0.0 

PINK(2) = 0.0 

DO 56 +L = 1,NAIRFO 

YMID = PIVOT * SINALF(L) + (L-1)*YSHIFT 
XLE = -PIVOT*COSALF(L)+( L-1)*XSHIFT 
aL = XLE 

ab = 0.0 

NPHI = 10 * NLOWER 

PDUM(L) = 0.0 

DO 30. I = 1,NPHI 

FRACT = FLOAT(I)/FLOAT(NPHI) 

XLP = -100.0 *TRY * (1.0 - COS(0. 5*PI*FRACT) )+XLE 
IF (I .EQ. 1) XLP = -0.000197+XLE 

XLP = -10.0 * (1.0 - COS(0. 5*PI*FRACT) ) 
DELX = XL - XLP 

XMID = 0. 5*(XL+XLP) 

XMID = 0. 5*(XL+XLP)*COSALF(L) 

YMID = 0.5*(XL+XLP)*SINALF(L) 

XL = XLP 

VELX = UGUST 


ADD CONTRIBUTION OF J-TH PANEL 


DO 40 MM = 1,NAIRFO 

LJ =(MM-1)**NODTOT 
KJ =(MM-1)*NP1 

DO 20 J = 1,NP1 

IF (J. EQ. NES GOTO} 24 


DXJ = XMID - X(J+KJ) 

DXJP = XMID - X(J+KJ+1) 

DYJ = YMID - Y(J+KJ) 

DYJP = YMID - Y(J+KJ+1) 

FLOG © = . 5**ALOG( (DX JP**DXJP+DYJP*DYJP) /(DXJ*DXJ+DYJI*DYJ) ) 
FTAN = ATAN2(DYJP*DXJ-DXJP*DYJ , DXJP*DXJ+DYJP*DYJ) 
CALMTJ = -COSTHE(J+KJ) 

SALMTJ = SINTHE(J+KJ) 

CALMTJ = -COSALF(L)*COSTHE(J+KJ) - SINALF(L)*SINTHE( J+KJ) 
SALMTJ = -SINALF(L)*COSTHE(J+KJ) + COSALF(L)*SINTHE( J+KJ) 
APY = PI2INV*(FTAN*CALMTJ + FLOG*SALMTJ) 

BPY = PI2INV*(FLOG*CALMTJ - FTAN*SALMTJ) 

VELX = VELX - DPROD(BPY,QK(J+LJ)) + DPROD(GAMK(MM) ,APY) 
GO TO 20 

DXJ = XMID - X(J+KJ) 

DXJP = XMID - X(2*NP1+MM) 

DYJ = YMID - Y(J+KJ) 

DYJP § = YMID - Y(2*NP1+MM) 

FLOG = .5*ALOG((DXJP*DXJP+DYJP*DYJP) /(DXJ*DXJ+DYJ*DYJ) ) 
FTAN = ATAN2(DYJP*DXJ-DXJP*DYJ, DXJP*DXJ+DYJP*DYJ) 
CALMTJ = -COSALF(L)*COSTHE(J+KJ) - SINALF(L)*SINTHE(J+KJ) 
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20 
40 


60 
70 
50 


gee 
30 


56 


230 


250 


240 


SALMTJ 
APY 
DUMMY 
VELA 
CONT TINUE 
CONTINUE 


- SINALF(L)*COSTHE( J+KJ) + COSALF(L)*SINTHE( J+KJ) 
PI2INV*(FTAN*CALMTJ + FLOG**SALMTJ) 

SS(MM)**( GAMMA( MM) -GAMK(MM) ) /WAKE( MM) 

VELX + DPROD( APY, DUMMY) 


ADD CORE VORTEX CONTRIBUTION 
remem EQ. 1) GO TO 50 


MM1 = 2-1 

DO 70 II = 1,NAIRFO 

KN = tT reat 

DOwGO0 oN = 1,MMI 

DX = XMID - XC(II,N) 

DY = YD ane (le NN) 

DIST = SQRT( DX*DX+DY*DY) 

COSTHN = DX/DIST 

SINTHN = DY/DIST 

SALMTN = -SINALF(L)*COSTHN + COSALF(L)*SINTHN 
Cer = -PI2ZINV*SALMIN/DIST 

VELX = VELX + DPROD(CPT,CV(II,N)) 
CONTINUE 

CONTINUE 


PDUM(L) = PDUM(L) + VELX * DBLE(DELX) 
FORMAT (215,6E14. 6,2F10. 6) 

CONTINUE 

DP(L) = PDUM(L)-PINK(L) 

CONTINUE 

PINK(1) = PDUM(1) 

PINK(2) = PDUM(2) 


COMPUTATION OF THE VELOCITY POTENTIAL FOR MIDPOINT OF EACH PANEL 


DO 240 L = 1,NAIRFO 
LI = (L-1)*NODTOT 
dele = -PINK(L) 


BEGIN WITH LOWER SURFACE 


DO 230 I = NLOWER,1,-1 
PHC = FHP*rHIkCItLi) 
PrexCI+L1) = 0. 5*( PHP+PHC) 
PHP = PHC 

PHITEL(L) = PHC 


RESET FOR UPPER SURFACE 


PHP = -PINK(L) 

DOs 250 I = NLOWER+1 ,NODTOT 
PHC = PRPtPHIK(CI+L1) 
PHIK(I+LI) = 0. 5*(PHP+PHC) 
PEP = PHC 

PHITEU(L) = PHC 


CONTINUE 
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C COMPUTE CP AT MID POINT OF I-TH PANEL 
C 

DO 295 L = 1,NAIRFO 

LI = (L-1)*NODTOT 

KI = (L-1)*NP1 

SUMC(L) = SUMC(L)/SS(L) 

DO 290 I = 1,NODTOT 

XIMID = .5*(XICI+K1) + XICIFKE 

XMID = . S(XCIRD) + XCI+tKI+T )) 

YMID =", SCG) eC th iD) 


CP(I+LI)= CP(I+LI) - 2.*(PHIK(I+LI)-PHI(I+LI))/DT 
WRITE (6,1050) I+LI,XIMID,XMID,YMID,QK(I+LI) ,GAMK(L) ,CP(I+LI), 
+ UE(I+LI),VN(I+L1) , PHIK(I+LI) ,PHI(I+LI) ,SUMC(L) 


290 CONTINUE 


WRITE (6,235) PINK(L) 


295 CONTINUE 


235 FORMAT (1X,'PHI AT LEADING EDGE =' ,F10.6,/) 


1000 FORMAT( / ,4%,-9 GX Ni eee 
CP(J)' 
INTGAMK’ ,/) 


+'GAMMA' ,4X 


+'PHI', 6x,! 


‘YC 6X MOGI. 6X, 


7X, ' VJ)" 6X, 'VNCJ)' ,5X, 'PHIK’ ,6X, 


1050 FORMAT(1I5,12F10.5) 


1200 FORMAT(1X,' LENGTH OF LEADING EDGE INTEGRATION IN CHORD =' ,F10. 6/) 


RETURN 
END 


cCcccccccccccccccccccccccccccccccccccccccccccccccccccccccCccccccccccccccc 


C2 Core G) C2 


SUBROUTINE CORVOR (SINALF ,COSALF) 
COMPUTE THE LOCAL VELOCITIES OF CORE VORTICES 


OaQaan 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 


GEVX(I Nee 
CCVY(I,N) .... 
XC(1,N) 
YC(1,N) 


GAMY(T) 


GBMY(1) 


AMY( 1) 


Bier) 


SUMAMY( 1 tegen. 


OOOO OO 2 OO © (Ore Gm 11 0) 0) 02: .) Cle oD 


X-VELO*OF THE I-TH AIRFOIL, N-TH CORE VORTEX Wii 
RESPECT TO THE CURRENT FROZEN FRAME OF REFERENCE 
Y-VELO OF THE I-TH AIRFOIL, N-TH CORE VORTEX WITH 
RESPECT TO THE CURRENT FROZEN FRAME OF REFERENCE 
X-COORD OF THE LOCATION OF THE I-TH AIRFOIL, N-TH 


CORE VORTEX W.R.T. THE GLOBAL FRAME OF REFERENCE. 


Y-COORD OF THE LOCATION OF THE I-TH AIRFOIL, N-TH 


CORE VORTEX W.R.T. THE GLOBAL FRAME OF REFERENCE. 


GLOBAL Y-VELO AT A LOCATION OF THE I-TH AIRFOIL 
CORE VORTEX DUE TO A SOURCE DIST OF UNIT STRENGTH 
ON ONE PANEL 

GLOBAL Y-VELO AT A LOCATION OF THE I-TH AIRFOIL 
CORE VORTEX DUE TO A VORTEX DIST OF UNIT STRENGTH 
ON ONE PANEL 

LOCAL Y-VELO AT A LOCATION OF THE I-TH AIRFOIL 
CORE VORTEX DUE TO A SOURCE DIST OF UNIT STRENGTH 
ON ONE PANEL 

LOCAL Y-VELO AT A LOCATION OF THE I-TH AIRFOIL 
CORE VORTEX DUE TO A VORTEX DIST OF UNIT STRENGTH 
ON ONE PANEL 

LOCAL Y-VELO AT A LOCATION OF THE I-TH AIRFOIL 
CORE VORTEX DUE TO A SOURCE DIST OF UNIT STRENGTH 
ON ALL PANELS 
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SUMBMY(I) .... LOCAL Y-VELO AT A LOCATION OF THE I-TH AIRFOIL 
CORE VORTEX DUE TO A VORTEX DIST OF UNIT STRENGTH 
ON ALL PANELS 


GCMX -.-» GLOBAL X-VELO AT A LOCATION OF A CORE VORTEX DUE 
TO ANOTHER POINT VORTEX OF UNIT STRENGTH 

GCMY .... GLOBAL Y-VELO AT A LOCATION OF A CORE VORTEX DUE 
TO ANOTHER POINT VORTEX OF UNIT STRENGTH 

CMX ..-. LOCAL X-VELO AT A LOCATION OF A CORE VORTEX DUE 
TO ANOTHER POINT VORTEX OF UNIT STRENGTH 

CMY .... LOCAL Y-VELO AT A LOCATION OF A CORE VORTEX DUE 


TO ANOTHER POINT VORTEX OF UNIT STRENGTH 


SUBROUTINE CORVOR 
COMMON /BOD/ IFLAG,NLOWER,NUPPER,NODTOT,X( 202) ,Y( 202), 
COSThEGCZ0 Pol NGheG 201). S502 )eNF] ,NP2Z,NPS,NP4, 
NP5 ,XSHIFT, YSHIFT ,NAIRFO,XI(202) ,YI(202), 
COSTHL( 201) , SINTHL( 201) 
COMMON /SING/ Q( 200) ,GAMMA( 2) ,QK( 200) ,GAMK( 2) 
COMMON /WAK/ VYW(2),VXW(2) ,WAKE(2) , DT 
Geren /CORV/ CVC2Z, 200) xC( 2, 200), YC( 2,200) ,M,TD,CCVX(2, 200), 
ar CON YC2e200)., XCIC2, 200), YCIC 2,200) 
COMMON /POT/ PHI(200) ,PHIK(200) 
COMMON /GUST/ UG(200),VG( 200) ,XGF,UGUST , VGUST 
COMMON /NUM/ PI,PI2INV 
COMMON /GEOM/ SINALF(2),COSALF(2) ,OMEGA( 2) ,UX(2),UY(2),PIVOT, 
~ APR PRM 
DIMENSION AMY( 2) , BMY(2) ,SUMAMY(2) ,SUMBMY(2), 
+VRC2), VYC2Z) pGAMYC2) ,GBNYC2) 
ir (MEQ, 1), GOTO 46 
MM1 =M- 1 


+++ 


ACCOUNT FOR THE FREE STREAM INCLUDING GUST EFFECTS 


UGC = 0.0 

VGC = 0.0 

DO 15 I = 1,NAIRFO 

KN = (I-1)*Mil 

DO 10) =6N = 1,Mfl 

XG = XC(CI,N)*COSALF(I) + YCCI,N)*SINALFCT) 
me=@CxG .GI. XGF) GO TO 5 

UGC = UGUST 

VGC = VGUST 

CONTINUE 

VY(I) = (1.+UGC)*SINALF(I)+VGC*COSALF(T) 
VXCI) = (1. +UGC)*COSALF(1)-VGC*SINALF(T) 
AMID = XCCI,N) 

YMID = YCCI,N) 


CALCULATE THE INFLUENCE COEFFICIENT DUE TO THE AIRFOILS 


DO 25 L = 1,NAIRFO 


LJ =(L-1)*NODTOT 
KJ =(L-1)*NP1 
SUMAMY(I) = 0.0 


SUMBMY( I) 0.0 
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Wl 


2 


20 


Zo 


50 
52 


10 


DO 20m Saat Tel 

DXJ = XMID - X(J+KJ) 
DYJ = YMID - Y(J+KJ) 
IF (J ,EQ, NPP) Geno 11. 


DXJP = XMID - X(J+KJ+1) 

DYJP = YMID - Y(J+KJ+1) 

GOunOr12 

DXJP = XMID - X(2*NP1+L) 

DYJP = YMID - Y(2*NP1+L) 

FLOG = .5*ALOG((DXJP*DXJP+DYJP*DY JP) /(DXJ*DXJ+DYJ*DYJ) ) 
FTAN = ATAN2(DYJP*DXJ-DXJP*DYJ, DXJP*DXJ+DYJP*DYJ) 
GAMY(I) = PI2INV*(FTAN*COSTHE(J+KJ) - FLOG*SINTHE( J+KJ)) 
GBMY(1I) = PI2INV*(FLOG*COSTHE( J+KJ) + FTAN*SINTHE( J+KJ) ) 


AMY(I) = GAMYCI)*COSALF( 1)-GBMY(C1I)*SINALF(I) 
BMY(1) = GAMY(1I)*SINALF(I)+GBMY(1I)*COSALF(T) 
IF (J. EQ: NPI eGeieeze 


SUMAMY(I) = SUMAMY(I) + AMY(TI) 
SUMBMY(I) = SUMBMY(I) + BMY(I) 
VY(I) = VY¥(I) + AMY(1)*QK(J+LJ) 
VX(I) = VX(I) - BMY(I)*QK(J+LJ) 
CONTINUE 


ADD THE CONTRIBUTION DUE TO THE WAKE ELEMENTS 


VY(I) = VY(I) + SUMBMY(I)*GAMK(L) + SS(CL)*BMY(1)* 
+ 


- 


(GAMMA( L) -GAMK(L) ) /WAKE(L) 

VX(I) = VX(I) + SUMAMY(I)*GAMK(L) + SS(L)*AMY(1)* 
(GAMMA(L) -GAMK(L) ) /WAKE(L) 

CONTINUE 


CALCULATE INFLUENCE COEFFICIENT DUE TO THE WAKE 


DO 35 L = 1,NAIRFO 

KMC = (L-1)*MM1 

DO 30 MC = 1,MM1 

IF (CL. EQ.1) -AND. (MC. EQN) Gono meo 


DX = XMID - XC(L,MC) 

DY = YMID - YC(L,MC) 

DIST2 = DX*DX+DY*DY 

GCMY = -PI2INV*DX/DIST2 

GCMX = +PI2INV*DY/DIST2 

CMY § = GCMY**COSALF(1)+GCMX*SINALF(I) 
CMX = -GCMY*SINALF( I )+GCMX*COSALF(I) 


ADD CORE VORTEX CONTRIBUTION 


VYCI) = VYC1) 4 GM Gy Cre) 
VXCI) = VXCI) + CMX*CVCL,MC) 


CONTINUE 


CONVECTION VELOCITY OF CORE VORTICES AT NEXT TIME STEP 


CCVX(I,N) = VX(I) 
CCVY(I,N) = VY(I) 
CONTINUE 


eZ 


15 CONTINUE 
40 CONTINUE 


RETURN 
END 
CGwerwoccec lu COCCCECCECCCCCCCCCCCCCCCECGGCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C C 
G SUBROUTINE NEWPOSN GC 
CG (e 
G TRANSFORM THE RESPECTIVE FROZEN LOCAL COORDINATES AND PANEL C 
C ANGLES TO THE GLOBAL FRAMES OF REFERENCE G 
& C 
ieee eececcCCOCCGCGCCCCECCGCCCECECCCCCCCCCCCCCCCCCCCCCECGCCCCCCCCCCCCCCCC 
C 
C XSHIFT .... GLOBAL X LOCATION OF THE 2-ND AIRFOIL PIVOT POSITION 
C YSHIFT .... GLOBAL Y LOCATION OF THE 2-ND AIRFOIL PIVOT POSITION 
G ITYPE FLAG FOR DEALING WITH THE VARIOUS COMPONENTS 
C = 0: FOR TRANSFORMING AIRFOIL COORDINATES 
C = 1 : FOR TRANSFORMING WAKE ELEMENTS 
G = 2 : FOR TRANSFORMING MOST RECENT CORE VORTEX SHED 
G = 3 : FOR TRANSFORMING ALL PREVIOUS CORE VORTICES 
C XPRM X% GLOBAL RELATIVE MOVEMENT OF THE PIVOTS POSITION 
¢ YRRM Y GLOBAL RELATIVE MOVEMENT OF THE PIVOTS POSITION 
C 
C 
SUBROUTINE NEWPOSC ITYPE) 
oo /BOD/ IFLAG,NLOWER ,NUPPER,NODTOT ,X( 202) ,Y( 202), 
COSTHE( 201), SINTHE( 201), $S(2), NP1,NP2,NP3,NP4, 
i NP5 ,XSHIFT,YSHIFT,NAIRFO X1( 202), YI(202), 
15 COSTHL( 201), SINTHL( 201) 
C@rmen /CORV/, CV(C2,200) ,XC( 2,200), YC(2;,200) ,M,TD,CVVxAC 2,200), 
+ Gree 200) xC1C2,200 ), YETEZ. 200 ) 
COMMON /GEOM/ SINALF(2) ,COSALF(2) ,OMEGA( 2) ,UX(2),UY(2), PIVOT, 
+ APRine eRe 
G 
G CHECK FLAG TO ASCERTAIN TYPE OF TRANSFORMATION 
C 
Pee(sliree . BO. 1 ) GOTO 30 
[Me ( ITYPE .—E©0. 2 ) GOTO 40 
TEC Were ont. 3 ) GOTO 45 
G 
c AIRFOIL COORDINATES AND LOCAL ANGLES TRANSFORMATION 
C 
DO 10 I = 1,NP1 
X( 1) = (XI(1)-PIVOT)*COSALF( 1)+YICI)*SINALF( 1) 
¥CT) = (XI(1)-PIVOT)*( -SINALF(1))+YIC1I)*COSALF( 1) 
10 CONTINUE 
DO 20 I = NP1+1,2*NP1 
X( 1) = (XI(1)-PIVOT)*COSALF(2 )+YIC 1I)*SINALF( 2 )+ 
- XSHIFT + XPRM 
rm 1) = (XI(1I)-PIVOT)*( -SINALF( 2) )+YIC1)*COSALF( 2 )+ 
+ YSHIFT + YPRM 
20 CONTINUE 
C 
G TRANSFORM LOCAL ANGLES 
C 
Dore 22eG = 1,2 
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ep ee ke 


kc) 2 


i I Se le 


fF) ©), C) 


P22 


30 


ae 


ae 


225 


40 


ze 


ak 


45 


4. 


+. 
46 
50 


K = (L-1)*(NODTOT+1) 
DO 222 J = 1,NODTOT 

DX = K(k ene IK ) 
DY = Y(J+1+K) - Y(J+K) 
DIST = SQRT(DX*DX +DY*DY) 


SINTHE( J+K) = DY/DIST 
COSTHE( J+K) = DX/DIST 
COTO SC 


WAKE ELEMENT TRANSFORMATION 


I = NP2 

AC®) = (XI(I)-PIVOT)*COSALF(1)+YIC I)*SINALF( 1) 

SIL, = (XI(1I)-PIVOT)*( -SINALF( 1) )+YICI)*COSALF( 1) 

I = NP2 + 1 

X(T) = (XI(1I)-PIVOT)*COSALF(2)+YICI)*SINALF( 2)+ 
ASHIFI + APR 

Al), = (XI(1I)-PIVOT)*(-SINALF( 2) )+YIC1I)*COSALF(2)+ 


YSHIPIT = VER 
WAKE ELEMENT ANGLES TRANSFORMATION 


DO 223 L = 1,2 
K (L-1)*(NODTOT+1) 


J = NP1 

DX = X(NP2+L-1) - X(NP1*L) 
DY = YCNP2+L-1) = YCNP ia) 
DIST = SQRT(DX*DX +DY**DY) 


SINTHEC IK )a—sDy isn 
COSTHE( J+K) = DX/DIST 
GOTO 50 


MOST RECENT CORE VORTEX SHED TRANSFORMATION 


XC(1,M) = (XCI(1,M)-PIVOT)*COSALF(1)+YCI(1,M)*SINALF(1) 

YC(1,M) = (XCI(1,M)-PIVOT)*( -SINALF( 1) )+YCI(1,M)*COSALF( 1) 

XC(2,M) = (XCI(2,M)-PIVOT)**COSALF( 2)+YCI(2,M)*SINALF(2)+ 
XSHIFT + XPRM 

YC(2,M) = (XCI(2,M)-PIVOT)*( -SINALF( 2) )+YCI(2,M)*COSALF( 2)+ 
YSHIFT + YPRM 

GOTO 50 


ALL PREVIOUS CORE VORTICES TRANFORMATION 


DO 46 I = 1,M-1 

XC(1,I) = (XCI(C1,1)-PIVOT)*COSALF( 1)+YCI(1,1)*SINALF(1) 

YC(1,I) = (XCI(1,1)-PIVOT)*( -SINALF( 1) )+YCI(1,1)*COSALF(1) 

XC(2,1) = (XCI(2,1)-PIVOT)*COSALF(2)+YCI(2,1)*SINALF(2)+ 
XSHIFT + XPRM 

YC(2,I) = (XCI(2,1I)-PIVOT)*(-SINALF( 2) )+YCI(2,1)*COSALF(2)+ 
YSHIFT + YPRM 

CONTINUE 

CONTINUE 

RETURN 

END 
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APPENDIX B. EXAMPLE INPUT DATA FOR PROGRAM USPOTF2 
5 


TVET Pee Pele tee te Be To Rete RI RICE Ra aR: aaa RRR elcktehdockteleseteletelots 


AIRFOIL TYPE : 8.4% THICKNESS SYMMETRICAL MISES AIRFOIL 
AIRFOIL COORDINATES ARE USER INPUT DATA (IFLAG = 1) 


NLOWER = 25 , NUPPER = 25 
tevetetede dedeteselevetee se tetek cetek kissd kcckhkkhdKcccKnkRecsickhdicntkcikhhkhcccdek 


01 /5S aS 
2 0.0 2.0 

1.000000 0.994858 0.980866 0.958884 0.929536 0.893455 
Ce on06,) OFGCSG IS = OF 751753 “0.695946 0.637271 0.576620 
0.514918 0.453098 0.392084 0.332794 0.276105 0.222865 
Ome col §O, 129819 O2091393 0.059146. 0.033560 0.015010 
0.003767 0.000000 0.003767 0.015008 0.033560 0.059146 
Seems O. 129619 OF17356)] 0.222865 0.276105 0.332791 
Oreee0e2) 0.453095 025514915 0.576617 0.637266 0.695946 
Ceyen7s0 0.603815 05851508 0.893455 0.929536 0.958884 
0.980866 0.994858 1.000000 

Uaeeeu00 -O. 000782 -0. 002754 =-0.00572Z1 -0.009351 -0.013459 
semouwoo7 -0,022285 -0O.026618 -0.030671 ~-0.034289 -0.037341 
pOmong7 12 =0.041315 -0,062083 -0.041979 -0.040979 -0.039096 
mUmeoteoUN-0.032820 -0. 0286555 -0.023651 =0..018220 -0.012379 
=Onevoz>7, 0.000000 0:,006259 0.012379 “0.018220 0.023651 
OmOzageoo 0.032820 0.056360 0.039096 0.040979 0.041979 
0.042083 0.041314 0.039712 0.037341 0.034289 0.030671 
Umeecors) 0.022755 OF01755/ 0.013459 0.009351 0.005721 
0.002784 0.000782 0.000000 

1.000000 0.994858 0.980866 0.958884 0.929536 0.893455 
Omowees 0.505815 0.759753 0.695948 0.637271 0.570620 
0.514918 0.453098 0.392084 0.332794 0.276105 0.222865 
tleecols 0.129819 0.091393 0.059146. 0.033560 0:015010 
Dyees767 0, 000000 0.003767 0.015008 0.033560 0.059146 
Gees §©.129819 @©2173661 6. 222865 0.276105 0. 332791 
Ome 2062) 90.653095 05574915 0.576617 0.637266 0.695946 
Orel 50720, 803815 0.851306 0.893455 0.929536 0.958884 
0.980866 0.994858 1.000000 

0.000000 -0.000782 -0.002784 -0.005721 -0.009351 -0. 013459 
=(-0lveee -0.022285 ~-0.026618 -0:030671 -0.034289 -0.037341 
-0.039712 -0.041314 -0.042083 -0.041979 -0.040979 -0. 039096 
~Om@ecmce) ~On037520 -O08028555 -0.023651 -0.018220 -0.0123/9 
—0gees) OG, 000000 OP006259 OF70T2379""0.018220 0.023651 
OmUgesses O1022820 OF086360 07089096". 040979 0.041979 
O7042065 07060514 0,039712 0.037341 0.034289" 0.030671 
GyOeztor: O2022285 0.017837 0.013459 0.009351 0.005721 
0.002784 0.000782 0.000000 

0.0 0.0 ~45.80000 0.00 0.0 . 0000 0 

0.0 0.00 0.0 0. 00 0.0 . 000 0.0 

5. 0000 0.065 0.001 OePOC feos 5 O9omee. 921370" 1. 2FO290 


135 


APPENDIX C. EXAMPLE OUTPUT DATA FROM PROGRANI 
USPOTF2 


DATA READ FROM FILE CODE 1 


5 

FIRB JOROBBORU BUBB ERI EB OUBO URE RU GE IEIOIEH IIE 
AIRFOIL TYPE : 8.4% THICKNESS SYMMETRICAL MISES AIRFOIL 
AIRFOIL COORDINATES ARE USER INPUT DATA (IFLAG = 1) 
NLOWER = 25 » NUPPER = 25 

FIBRO EUR ERO BUBB BBO BB BIBER BEBE BEBE 


oo oo o0o0o00 F 


1 25 
2 0.0 
-000000 
-851308 
-514918 
- 173861 


.003767 
-091393 
- 392082 
- 751750 
- 980866 
-000000 
-0. 
-039712 


017837 


.036360 
-006259 
-028555 
-0420835 
.026618 


.00278¢ 


.000000 
.851308 
-514918 
- 173861 
.003767 
-091393 
- 392082 


- 751750 
- 980866 
-000000 


-017837 
.039712 
.036360 
.006259 
.028555 
-042083 
.026618 
-002784% 
.000000 
-000000 
-650000 


25 
2.0 
0.994858 
0.803815 
0.453098 
0.129819 
0.000000 
0.129819 
0.453095 
0.803815 
0.994858 
-0.000782 
-0.022285 
-0.041314 
-0.032820 
0.000000 
0.032820 
0.041314 
0.022285 
0.000782 
0.994858 
0.803815 
0.453098 
0.129819 
0.000000 
0.129819 
0.453095 
0.803815 
0.994858 
-0.000782 
-0.022285 
-0.041314 
-0.032820 
0.000000 
0.032820 
0.041314 
0.022285 
0.000782 


0.980866 
0.751753 
0.392084 
0.091393 
0.003767 
0.173861 
0.519915 
0.851308 
1.000000 
-0.002784 
-0.026618 
-0.042083 
-0.028555 
0.006259 
0.036360 
0.039712 
0.017837 
0.000000 
0.980866 
0.751753 
0.392084 
0.091393 
0.003767 
0.173861 
0.514915 
0.851308 
1.000000 
-0.002784¢ 
-0.026618 
-0.042083 
-0.028555 
0.006259 
0.036360 
0.039712 
0.017837 
0.000000 


0.000000-45.800003 


0.000000 
0.025000 


0.000000 
0.001000 


0.958884 
0.695948 
0.332794 
0.059146 
0.015008 
0.222865 
0.576617 
0.893455 


-0.005721 
-0.030671 
-0.041979 
-0.023651 
0.012379 
0.039096 
0.037341 
0.013459 


0.958884 
0.695948 
0.332794 
0.059146 
0.015008 
0.222865 
0.576617 
0.893455 


-0.005721 
-0.030671 
-0.041979 
-0.023651 
0.012379 
0.039096 
0.037341 
0.013459 


0.000000 
0.000000 
0.000000 


0.929536 
0.637271 
0.276105 
0.033560 
0.033560 
0.276105 
0.637266 
0.929536 


-0.009351 
-0.034289 
-0.040979 
-0.018220 
0.018220 
0.040979 
0.034289 
0.009351 


0.929536 
0.637271 
0.276105 
0.033560 
0.033560 
0.276105 
0.637266 
0.929536 


-0.009351 
-0.034289 
-0.040979 
-0.018220 
0.018220 
0.040979 
0.034289 
0.009351 


0.000000 


0.000000 
4.855698 
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0.893455 
0.576620 
0.222865 
0.015010 
0.059146 
0.332791 
0.695946 
0.958884 


-0.013459 
-0.037341 
-0.039096 
-0.012379 
0.023651 
0.041979 
0.030671 
0.005721 


0.893455 
0.576620 
0.222865 
0.015010 
0.059146 
0.332791 
0.695946 
0.958884 


-0.013459 
-0.037341 
-0.039096 
-0.012379 
0.023651 
0.041979 
0.030671 
0.005721 


0.000000 0.000000 0 
0.000000 0.000000 0.000000 
-0.921370 1.216290 0O 


lec 


826000°O0 028600°O O000000°0 »S92280'°T- bZ2TZT°O- O£6000°0 S£E8Sb0°0- 9256850°0- 692565'°0 2T 
826000°0 S086T0°O O00000°0 ThzbZ0°T- 8025ST°O- 0£6000°0 »908S50'°0- ST8S"50'°0- S66909°0 IT 
826000°0 [S5T6T0°O O00000°O0 [22S590'T- 606b<T*O- 0£6000°0 Z82690°0- O8b2Z20'°0- 609999°0 OT 
826000°O0 <£08220°O TO00000°O Z£0SSO*T- 2OTSTT’0- 0£6000°0 90080'0- %59820'0- OS8EZZ'0 6 
826000°0 2699520°0 TO00000'O0 %22<2b0°T- 026880°0- 0£6000°0 656680'0- IS>62Z0'0- 4842/2°0 8 
826000°0 692220°0 TOO0000°O0 229620°T- S£2090°0- O0£6000°0 O022660°0- T900Z0°0- 19S5/28'0 Z 
826000°0 020620°O0 TO00000°O0 T89ETO'T- 055220°0- O£6000°0 T[S6ZOT*O- Sb9STO’O- Tes2Z8°0 9 
826000°0 »94620°0 TO00000°O0 965666'°0- 8ZZ40TO°O O£6000°0 S829STT’°O- SObTIO'O- S65TT6°0 SS 
826000°O SZT620°O0 000000°0 Ob0TZ6°0- T80ZS50°0 O£6000°0 TI2Z22T’0- 9£59/00'0- OTZ5%6'°0 » 
826000°O0 T62820°O TOO0000'O- £55056°0- BSbSSTI*°O O£6000°0 68b6L2T°0- 292500°0- G28696'0 ¢£ 
826000°0 660220°O0 TO0000'0- 906268°0- 299b6T°O O£6000°0 80962T’*0- £8ZT00°0- 2982986'°0 2 
826000°O 126520°O TO00000°0- 605628°0- ST6TT£°O O£6000°0 09292T*O- T6éZ000°0- 626266'°0 T 
VWAVSINI IHd (CT ONA (TIA (f ddd VWWYS (FD CTIA CTX c 

ic Y3SGWNN TIOsYTV 

000000°0 = (2 JVHdIV LY NOIINIOS MOIS AGV3IS 

000000°0 = €T JVHd1V LY NOTiNIOS MO14 AGV31S 

0°2 = J4SIHSA 0°0O = L4IIHSX 

66S8T0'2 = HLON3I] YILSWIN3d (2 )IIOINIV 

66S58T0°2 = HISN317 SILSWIN3d (T J TIOASUTY 


Ni 
0 


SVIOSUTV 40 ON 1VLOL 


826000°0 
826000°0 
826000°0 
826000°0 
826000°0 
826000 °0 
826000°0 
826000°0 
826000°0 
826000°0 
826000°0 
826000°0 
826000°0 
826000°0 
826000°0 
826000°0 
826000°0 
826000°0 
826000°0 
826000°0 
826000°0 
826000°0 
826000°0 
826000°0 
826000°0 
826000°0 
826000°0 
826000°0 
826000°0 
826000°0 
826000 °0 
826000°0 
8260000 
826000°0 
826000°0 
826000°0 
826000°0 
826000°0 


<08220°0 
816820°0 
TST020°0O 
9TOTZ0°O 
0Z22T20°O 
9$22020°0 
5S5620°0 
$82220°0 
$82520°0 
£05020°0 
STO9TO°O 
2260T0°0 
562900°0 
485000 °0- 
T2Z900°0- 
6T62T0°O- 
6206T0°0- 
9$06%20°0- 
z<T$020°0- 
926920 °0- 
628620 °O- 
869250 °0- 
105960 ° 0- 
Z£T9860°0- 
61T96%0°0- 
0196%0°0- 
469850°0- 
2259%0°0- 
929250 '0- 
$28620°0- 
696520 °0- 
T29020°0- 
660520°0- 
TZ26T0°O- 
8HzzZTO° O- 
62200°0- 
2TzT00°0- 
5S5500°0 


100000 °0- 
T00000°0- 
100000 °0- 
T00000°0- 
T00000°0 
100000°0 
T00000°0 
100000 °0 
T00000°0 
000000°0 
000000°0 
000000°0 
000000°0 
000000°0 
000000°0 
000000°0 
000000°0 
000000°0 
100000°0 
000000°0 
1T00000°0- 
~00000° 0- 
200000 °0- 
200000° 0- 
200000°0- 
000000°0 
T00000°0- 
TO00000°*0- 
100000 °0O- 
T00000° 0- 
000000°0 
T00000°0 
000000°0 
000000°0 
000000°0 
000000°0 
000000°0 
000000°0 


80S628°0 
229268°0 
28606 *0 
9zZ8TZ6°0 
£69966 °0 
866H5TO'T 
So2Tz0 'T 
46090 °T 
020250°T 
04S290°T 
869920 °T 
9TZ580°T 
»99T60°T 
199260°T 
S6Z20T'T 
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eT2Ortt tf 
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g26TZ20°T- 
019260°T- 
Z£o<Z0T*T- 
6860TT*TI- 
S29 UCL L= 
650LT*T- 
29080T °T- 
599H0T °T- 
£5z00T *T- 
59TS60°T- 
T9T680°T- 


9T6TTZ°O 

98Tbé6T "0 

2SSHTt’O 

$2S550°0 

2£89800°0 

T22020°0- 
996290 °0- 
022260'°0- 
9622TT°O- 
29962T O- 
62£269T °0- 
8099LZT*0- 
Oz<ZT6T’°O- 
658602 °0- 
9ST9T2°0- 
ZT9S22°0- 
TZS9222°0- 
2S62Z22 °0- 
£62622 °O- 
£92222 °0- 
622822 °0- 
996502 °0- 
92T65T *0- 
092820°0 

65S689°0 

T26989°0 

265920°0 

Oz T6%T *O- 
962502 °0- 
8T2922 °0- 
S62622°0- 
222S%2 °0- 
S6TZ22 °0- 
008222 °0- 
282022 ° O0- 
29L0T2 “0- 
S8z66T’0- 
TZ298T “O0- 


026000°0 
026000°0 
026000°0 
026000°0 
026000°0 
026000°0 
026000°0 
026000°0 
026000°0 
026000°0 
026000°0 
026000°0 
026000°0 
026000°0 
026000°0 
026000°0 
026000°0 
026000°0 
026000°0 
026000°0 
026000'°0 
026000°0 
026000°0 
026000°0 
026000°0 
026000°0 
026000°0 
026000°0 
026000°0 
0%6000°0 
026000°0 
026000°0 
026000°0 
026000°0 
026000°0 
026000°0 
026000°0 
026000°0 


282660° 


06982T °0- 
SOTT>T*O- 
VzT8zt'O- 
266TET*O- 
8Lbb2T'O- 
26ZSTT° O- 
9Z090T °0- 
29S960°0- 
T9980 °0- 
£29220°0- 
246650 °0- 
855950 °0- 
TOT220°0- 
2529T0°O- 
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